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ABSTRACT 

We describe the GALFORM semi-analytic model for calculating the formation and 
evolution of galaxies in hierarchical clustering cosm ologie s. It improves upon, and ex- 
tends, the earlier scheme developed by Cole et al. ( 1994 ). The model employs a new 
Monte-Carlo algorithm to follow the merging evolution of dark matter halos with ar- 
bitrary mass resolution. It incorporates realistic descriptions of the density profiles of 
dark matter halos and the gas they contain; it follows the chemical evolution of gas and 
stars, and the associated production of dust; and it includes a detailed calculation of 
the sizes of disks and spheroids. Wherever possible, our prescriptions for modelling in- 
dividual physical processes are based on results of numerical simulations. They require 
a number of adjustable parameters which we fix by reference to a small subset of local 
galaxy data. This results in a fully specified model of galaxy formation which can be 
tested against other data. We apply our methods to the ACDM cosmology (fio — 0.3, 
Ao — 0.7), and find good agreement with a wide range of properties of the local galaxy 
population: the B-band and K-band luminosity functions, the distribution of colours 
for the population as a whole, the ratio of ellipticals to spirals, the distribution of 
disk sizes, and the current cold gas content of disks. In spite of the overall success of 
the model, some interesting discrepancies remain: the colour-magnitude relation for 
ellipticals in clusters is significantly flatter than observed at bright magnitudes (al- 
though the scatter is about right), and the model predicts galaxy circular velocities, at 
a given luminosity, that are about 30% larger than is observed. It is unclear whether 
these discrepancies represent fundamental shortcomings of the model or whether they 
result from the various approximations and uncertainties inherent in the technique. 
Our more detailed methods do not change our earlier conclusion that just over half 
the stars in the universe are expected to have formed since z ^ 
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1 INTRODUCTION 



199e| ; |Ellis et al. 1996| ) and 2 ~ 3-4 ( ^teidel et al. 1999| ); the 



The past few years have been a remarkably rich period in 
observational studies of galaxy formation. Major advances 
have resulted from observations at many wavelengths, from 
the far ultraviolet to the sub-millimeter. Breakthroughs in- 
clude the discovery and measurement of the clustering of 
"Lyman-break" galaxies, a population of luminous, star- 
forming galaxies at red shifts 2 ~ 3 — 4 (Steidel et al. 1996 
Adelberger et al. 1998 ); estimates of the history of star for- 
mation and the attendant pro duction of metals, from z ~ 5 
to the present (Madau et al. 1996,1998); meas urements of 
the galaxy luminosity function at 2: ~ 0.5 — 1 (Lilly et al 



discovery of a population of bright sub-millimeter sources, 
some of which, at least, appear to be dusty, star-forming 
galaxies at z ^ 2 (Ivison et al. 1998). All of these and many 
other observations are beginning to sketch out an empirical 
picture of galaxy evolution. 

On their own, the data provide only a partial descrip- 
tion of specific stages of galaxy evolution. To develop a phys- 
ical understanding of the processes at work, and to relate 
observations to cosmological theory, requires detailed mod- 
elling that exploits our current understanding of astrophys- 
ical processes in their cosmological context. The theoretical 
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infrastructure required for this programme has b een in place 
for over a decade (e.g. Blumenthal et al. 1984; Pavis et al 



n its standard form, it assumes that galaxies grew 
out of primordial Gaussian density fluctuations generated 
during inflation and amplified by gravitational instability 
acting on cold dark matter, the dominant mass component 
of the Universe. Gas is initially mixed in with the dark mat- 
ter, and when dark matter halos collapse, the visible com- 
ponent of galaxies accumulates as stars condense out of gas 
that has cooled onto a disk. 

To construct a theory of galaxy formation that can be 
tested against observations requires combining the theory 
of the evolution of cosmological density perturbations with 
a description of various astrophysical processes such as the 
cooling of gas in halos, the formation of stars, the feedback 
effects on interstellar gas of energy released by young stars, 
the production of heavy elements, the evolution of stellar 
populations, the effects of dust, and the merging of galaxies. 
The most appropriate methodology is to carry out ab initio 
calculations that follow directly the development of primor- 
dial density fluctuations into luminous galaxies. Within the 
standard cosmological model, the initial conditions are very 
well defined. They are specified by the power spectrum of 
primordial density perturbations, whose shape is fixed by 
the cosmological parameters: the mean mass density, Oq, the 
mean baryon density, S7b, the cosmological constant, Aq, and 
the Hubble constant, Hq (which, throughout this paper, we 
express as Ho = 100/i kms~^Mpc~^.) 

The subsequent evolution of the dark matter and 
baryons is best calculated by Monte Carlo simulation. Two 
different approaches have been developed for this purpose. 
In the first, direct simulations, the gravitational and hydro- 
dynamical equations in the expanding universe are solved 
explicitly, using one or more of a variety of numerical tech- 
niques that have been specificall y developed for this purpose 
over the past twenty years (e 
al. 199^ 



Frenk et al. 1996 



Katz et al. 19921 [Eyrard et 



1999|; [Katz et al. 1996|; [Navarro 



Steiifmetz 1997|; pearce et al. 199E| ; p:hacker et al. 2000t 



Blanton^r^J~20o6[ ). In the second approach, now com- 



monly known as "semi-analytic modelling of galaxy forma- 
tion" ([White fc Rees 1978|; [White 



Kauffmann 



_ Frenk 1991| 

et al. 1 |)93| ; Cole et al. 1994), the evolution of the baryonic 
component is calculated using simple analytic models, while 
the evolution of the dark matter is calculated either directly, 
using N-body methods, or using a Monte-Carlo technique 
that follows the formation of dark matter halos by hierar- 
chical merging. It is this second approach that we discuss in 
this paper. 

The two modelling techniques have complementary 
strengths. The major advantage of direct simulations is that 
the dynamics of the cooling gas are calculated in full general- 
ity, without the need for simplifying assumptions. The main 
disadvantage is that even with the best codes and fastest 
computers available today, the attainable resolution is still 
some orders of magnitude below that required to resolve 
the formation and internal structure of individual galaxies 
in cosmological volumes. In addition, a phenomenological 
model, similar to that employed in semi-analytic modelling, 
is required to include star formation and feedback processes 
in the simulation. These processes are, in fact, much more 
difficult to treat and much more uncertain than the dynam- 
ics of the diffuse gas. 



Semi-analytic modelling does not suffer from resolu- 
tion limitations, particularly when Monte-Carlo methods are 
used to generate the halo merger histories. In this case, the 
resolution can be made arbitrarily high at a relatively small 
computational cost. The major disadvantage is the need for 
simplifying assumptions in the calculation of gas properties, 
such as spherical symmetry or a particular flow structure. 
It is encouraging that detailed comparisons between direct 
and s emi-an alytic simulati ons sh ow good agreement (Pearce 
et al. 1999 , Benson et al. 2000c ). An important advantage 
of semi-analytic modelling is its flexibility. This allows the 
effects of varying assumptions or parameter choices to be 
readily investigated and makes it possible to calculate a wide 
range of observable galaxy properties, such as luminosities 
in any waveband, sizes, mass-to-light ratios, bulge-to-disk 
ratios, circular velocities, etc. 

Semi-anal ytic m odellin g has its roots in t he wo rk of 



White & Rees (1976^Cole (1991), Lacey & Silk (1991), and 



White & Frenk (1991) who laid out the overall philosophy 
and basic methodology of this approach. Throughout most 
of the 1990s, this technique was developed and promoted 
primarily by two collaborati ons, one cu rrently base d at M u- 
nich (e.g. Kauffmann et al. 



Kauffm ann, 
19984b, 



1993 



Nusser & Steinmetz 



1994[ Kauffmann 
Mo 



199c, Kauffmann et a l. 1999a ) 



al. 



Du rham (e. g. Co le et al. 1994, Heyl e t al. 1995, Baugh et 



1997 



1995a ,b 



Mao & White 
and the other at 



19964 b, |l998| Benson et al. |2000 a; see also Lacey et al 



199J ). In the past two years, several other groups have begun 



to apply this technique to study vario us asp ects of galaxy 
for matio n (e.g. Avila-Reese & Firmani 199q, Guiderdoni et 



al. [19981, W u, Fab ian & Nulsen [1998 

199£, Somerville & Primack 



Peacock 



van Kampen, Jimenez 
This body 



199£) 



of work has demonstrated the usefulness of semi-analytic 
modelling as a means for fleshing out the observable conse- 
quences of current cosmological theories and for the inter- 
pretation of observational data, particularly at high redshift. 

A growing body of galaxy properties has been analysed 
using semi-analytic methods. Examples of noteworthy suc- 
cesses include the ability to reproduce the local fleld galaxy 
luminosity function, the slope and scatter of the TuUy-Fisher 
relation for spiral galaxies, and the counts and redshift dis- 
tributions of fai r it galaxies (e.g. see Kauffmann et al. 1993 



Cole et al. 1994; Kauffmann et al. 1994). Nevertheless, some 



important properties have remained obstinately difficult to 
reproduce, most notably the colour-magnitude relatio n for 
cluster ellipticals (but see Kauffmann & Chariot 1998a) and 
a simultaneous fit to the local luminosity function a nd th e 
zero-point of the TuUy- Fisher relation (e.g. Heyl et al. 1995). 

A wide variety of physical processes are involved in the 
formation of galaxies. Some of them, like star formation, are 
very poorly understood. Modelling galaxy formation there- 
fore inevitably requires making approximations and adopt- 
ing simplified descriptions of some of these processes. Most 
often, an incomplete understanding of a physical ingredi- 
ent is subsumed within a simple scaling law that contains 
free parameters. A remarkable facet of modern semi-analytic 
modelling is that a realistic picture of galaxy evolution can 
be formulated with a relatively small number of such pa- 
rameters, typically four or five in the simplest versions. A 
strategy that has proved useful is to fix the values of these 
parameters by trying to match a subset of local data (for 
example, the luminosity functions in two passbands or the 
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Figure 1. A schematic showing how different physical processes are combined to make predictions for the observable properties of 
galaxies, starting from initial conditions specified by the cosmology. The numbers in each box indicate the subsection of the paper in 
which our method for modelling that process is described. 



TuUy-Fisher relation). This leads to a completely specified 
model that has predictive power and may be used to calcu- 
late theoretical expectations for other local properties or for 
properties at high redshift. This approach h as me t with con- 
siderable success. For example, Cole et al. (1994) predicted 
that most of the stars in the universe formed at relatively 
low redshift [z ^ 1 for an Qo ~ I standard CDM cosmol- 
ogy) ; Kauffmann ( 1996| ) predicted a sharply declining num- 



ber of brigh t elliptical galaxies at hi gh re dshift; and Baugh 
et al. ( 1998 ) and Governato et al. (1995) predicted strong 
clustering for Lyman-break galaxies at redshift z ~ 3. 

In this paper we present a new semi-analytic mode l 
which builds upon the scheme described by Cole et al. (11994) 
which we us ed for a number of applications (Heyl et al.ll99q, 
Baugh et al. 1996a 3). Our new model differs from the earlier 
one primarily in its greater scope and richness, but also in 
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Figure 2. A schematic siiowing tlie observable galaxy properties predicted by the model. The numbers in each box indicate the subsection 
of the paper in which the comparison of model predictions with observations of that property are described. The predictions for galaxy 
clustering are described in separate papers, as indicated in the box. 



the manner in which certain key physical properties are cal- 
culated. These improvements are called for both by recent 
theoretical developments and, most importantly, by the in- 
crease in the quantity and quality of observational data. The 
main additions to our new code are the inclusion of chemical 
enrichment and dust processes, prescriptions for calculating 
the sizes of disks and spheroids, the use of more realistic den- 
sity profiles for dark matter halos and gas, and the ability to 
follow the mergers of halos with fine mass and time resolu- 
tion. It should be noted that despite all these improvements, 
when one adopts the same cosmological parameters and also 
the same galaxy formation parameters (e.g. for stellar feed- 
back), then the main predictions of the model, including the 
galaxy luminosity function, TuUy-Fisher relation and overall 
star format ion h istory, are practically identical to those in 
Cole et al. (1994). The one exception is the inclusion of dust 



extinction which makes t he gal axy colours somewhat redder 
than in the Cole et al. (1994) models. Thus, the changes 



to the model may be viewed as refinements that allow more 
properties of the galaxy population, such as galaxy sizes and 
metallicities, to be calculated. 

The main aim of this paper is to lay out the methods 
that we use in our new semi-analytic model and to compare 
results with a restricted set of observational data. This is 
a long paper containing a mixture of technical descriptions 
and results of more general interest. In the following, brief 



section, we present an overview, together with schematics 
illustrating how different parts of the model fit together. 
Non-specialist readers may wish to skip the more detailed 
passages of the paper on a first reading and we recommend 
how this might done in Section ^. 



2 OVERVIEW 

Our galaxy formation model is a synthesis of many tech- 
niques, each of which has been developed to treat particular 
aspects of the complex process of galaxy formation. Its back- 
bone is a Monte Carlo method for generating "merger trees" 
to describe the hierarchical growth of dark matter halos. The 
full range of properties and processes that we model within 
this framework are: 

(i) The gravitationally-driven formation and merging of 
dark matter halos. 

(ii) The density and angular momentum profiles of dark 
matter and shock heated gas within dense non-linear halos. 

(iii) The radiative cooling of gas and its collapse to form 
centrifugally supported disks. 

(iv) The scalelengths of disks based on angular momen- 
tum conservation and including the effect of the adiabatic 
contraction of the surrounding halo during the formation of 
the disk. 
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(v) Star formation in disks. 

(vi) Feedback, i.e. the regulation of the star formation 
rate resulting from the injection of supernova (SN) energy 
into the interstellar medium (ISM). 

(vii) Chemical enrichment of the ISM and hot halo gas 
and its influence on both the gas cooling rates and the prop- 
erties of the stellar populations that are formed. 

(viii) The frequency of galaxy mergers resulting from dy- 
namical friction operating on galaxies as they orbit within 
common dark matter halos. 

(ix) The formation of galactic spheroids, accompanied by 
bursts of star formation, during violent galaxy-galaxy merg- 
ers, and estimates of their effective radii. 

(x) Spectrophotometric evolution of the stellar popula- 
tions. 

(xi) The effect of dust extinction on galaxy luminosities 
and colours, and its dependence on galaxy inclination. 

(xii) The generation of emission lines from interstellar gas 
ionized by young stars. 

Our treatment of each of these processes is described in the 
foUowiiig sections. The model is summarized schematically 
in Fig. hi which also shows in which subsection of the paper 
each of the processes is described. 

The scheme we present is largely modular, and within 
each module one has the choice of selecting various options 
as well as certain parameter values. The options (for exam- 
ple, including or ignoring the baryonic mass of the galaxy 
when computing its rotation curve) are not degrees of free- 
dom within the model. Instead they allow us to vary the 
complexity of the description in order to gain physical un- 
derstanding. By switching processes on and off and changing 
certain assumptions we are able to determine which physi- 
cal processes are directly responsible for a particular galaxy 
property. 

The observable galaxy properties predicted by the 
model are shown schematically in Fig. H. This figure sep- 
arates the predicted quantities into two categories. Some 
aspects of the observational quantities in the first category 
are used as primary constraints on the model parameters. 
The boxes in the figure also indicate in which subsection of 
the paper the observational comparison for that quantity is 
presented, or in the case of galaxy clustering, the related pa- 
pers in which they are presented. Predictions for the redshift 
evolution of galaxy properties will be presented in future pa- 
pers. 

The layout of the paper is as follows: Section |^ presents 
techniques for generating merger trees to describe the grav- 
itational growth of dark matter halos and models for their 
internal structure. Section ^ describes how we calculate disk 
and spheroid formation, including star formation, feedback 
and chemical evolution, and the scalelengths of galactic disks 
and bulges. Section ^ presents our methods for calculat- 
ing the luminosities of stellar populations, and the effects 
of dust extinction. These three sections may be skipped on 
a first reading, simply noting the definition of the parame- 
ters that describe our star formation and feedback model. 



equations (4.4), (4.5), (4.14), and (4.15). An overview of 



the model and the strategy we adopt to constrain parame- 
ters and to obtain a well-specified, predictive model are pre- 
sented in Section ^, which should be ofgeneral interest. This 
procedure is implemented in Section M, where we illustrate 



the effects of varying each model parameter. The general 
reader may simply wish to study the figures in this section 
and read the summary in subsection 7.!;. In Section ^we test 
our fully specified model against further properties of the 
observed galaxy population. Again, the general reader may 
wish to study only the figures in this section. Finally, in Sec- 
tion 1^ we restate the general philosophy of our approach and 
discuss the strengths and weaknesses of the specific model 
that we have presented. This, and the final section present- 
ing a summary of our results, are both of general interest. 



3 FORMATION OF DARK MATTER HALOS 

Galaxies are assumed to form inside dark matter halos, and 
their subsequent evolution is controlled by the merging his- 
tories of the halos containing them. It is therefore essen- 
tial to have an accurate description of how dark halos form 
and evolve through hierarchical merging, and of the internal 
structure of these halos. These are both described in this 
section. 



3.1 Dark Matter Halo Merger Trees 

We use a new Monte-Carlo algorithm to generate merger 
trees that describe the formation paths of randomly selected 
dark matter halos. It is an impr ovement over th e "block 
model" that we used prev iously (Cole et al. 1994; Heyl et 



al. 1995; Baugh et al. 1996a ,b). The new algorithm is directly 
based on the analytic e xpress ion for halo merger rates de- 
rived by Lacey & Cole (1993). At each branch in the tree, 



a halo splits into two progenitors, but unlike in the "block 
model" , the mass ratio of the progenitors can take any value. 
Below, we briefly describe this new algorithm, and the way 
in which a population of merger trees is set up to provide a 
framework for modelling the processes of galaxy formation. 

It is possible to generate merger trees directly by follow- 
ing the evolution of dark matter halos in coUisionless cos- 
mological N- bodv simulations (e.g. Rou kema 



1997; KaufiF- 



mann et al. 1999a; van Kampen et al. 1999). Combining 
semi-analytic modelling and N-body simulations certainly 
provides a very powerful technique to investigate small 
scale galax y clustering (e. g. Ka uffmann, Nusser & Stein- 



metz 



1997, G overna to et al. 



1998 



Kauffmann et al. 



1999a ,b 



Benson et al. 2000a, b). However, extracting merger trees di- 
rectlty from N-body simulations carries the high price of a 
limited dynamic range in mass and much greater computa- 
tional complexity. Also, for many applications this appears 
to be unnecessary, as the properties of halo merger trees 
are u ncorrelated with environment (Lemson & Kauffmann 



199£), and the Monte Carlo merger trees agree well, statist! 
cally, with those extr a cted from N-body si r nulations ( [ Kauff- 



mann fc White 1995 ; Lacey & Cole 1994| jSomerville et al 
200C ; Lacey & Cole in preparation 



3.1.1 A New Algorithm 

Our starting point for generating a merger tree to describe 
the history of mergers experienced by an individual dark 
matter halo is equation (2.15) of Lacey & Cole (1993): 

1 {5cl — l5c2) 



/i2(Mi,M2)dMi = 



2^ (a2-<T2)3/2 
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X exp 



dMi.(3.1) 



This equat ion, derived from the extension of th e Pre ss & 
Schechter ( [L974 ) theory proposed by Bond et al. ( 1991 ) and 
Bower (1991), gives the fraction of mass, fi2{M\, M2)dMi, 
in halos of mass M2, at time t2, which at an earlier time, ii, 
was in halos of mass in the range Mi to Mi + dMi. Here, 
the quantities ai and (J2 are the linear theory rms density 
fluctuations in spheres of mass Mi and M2. The (5ci and 
5c2 are the critical thresholds on the linear overdensity for 
collapse at times ti and t2 respectively. (Specifically, a{M) 
and 5c{t) are the values extrapolated to z = according 
to linear theory.) For a critical density (SI = 1) universe, 
we adopt 5, — 1.686(1 + z), while for low Qo, open and 
flat, universes we adopt the appropriate expressions in the 
appen dices of Lacey & Cole (1993) and Eke, Cole & Frenk 
(1996) respecti vely. 

Equation ( |3.l[ ) can be used to estimate the recent 
merger histories of a set of halos w hich at time t2 have mass 
M2. Taking the limit of equation ( |3.l[ ) as ti — > t2, we ob- 
tain an expression for the average mass fraction of a halo 
of mass M2 which was in halos of mass Mi at the slightly 
earlier time ti, 



dti 



dMidti 



1 



dSci daf 



2^ (0-2 - a2)3/2 dti dMj 



dMidti. (3.2) 



Thus, the mean number of objects of mass Mi that a halo 
of mass M2 "fragments into" when one takes a step dti back 
in time is given by 



d/12 Afa 



(Ml < M2). 



(3.3) 



dN 

ImI 

This expression gives the average number of progenitors as 
a function of the fragment mass, Mi. It is this simple ex- 
pression that our algorithm uses to build a binary merger 
tree. 

The quantities that must be specified in order to define 
the merger tree are the density fluctuation power spectrum, 
which gives the function a(M), and the cosmological pa- 
rameters, f2o and Aq, which enter through the dependence 
of 5c{t) on the cosmological model. There is also one numer- 
ical parameter involved, the mass resolution, Mrcs. Having 
specified these parameters one can compute the two quan- 
tities, 



M2/2 



dN 
dWi 



dMi, 



(3.4) 



which is the mean number of fragments with masses. Mi , in 
the range M^cs < Mi < M2/2, and 



F = 



dN Ml 
dMi M2 



dMi, 



(3.5) 



which is the fraction of mass in fragments with mass be- 
low the resolution limit. Note that both these quantities are 
proportio nal t o the timestep dti , through the dependence in 
equation (3.3). 



Once these quantities have been defined, the algorithm 
to generate the merger trees is simple. First, choose a 
timestep, dti, such that P ^ 1, to ensure that multiple 
fragmentation is unlikely. Next, generate a random number, 
R, drawn uniformly from the interval to 1. If i? > P, then 



the main halo does not fragment at this timestep. However, 
the original mass is reduced to account for mass accreted in 
the form of halos with masses below the resolution limit, to 
produce a new halo of mass M2(l — F). If, however, R < P, 
then a random value of Mi in the range Mros < Mi < M2/2 
is generated, consistent with the distribution given by equa- 
tion (3.3), to produce two new halos with masses Mi and 
M2(l — F) — Ml. The same operation is repeated on each 
fragment at successive timesteps going back in time, and 
thus a merger tree is built up. 

The main advantage of this new al gorithm over the 
"block model" t hat we used previo usly ( Cole et al. 1994 



Heyl et al. 1995 



[Baugh et al. 199'6ab), and which has been 
used recently by Wu et al. ( 199^ , 200C ), is that there is no 
quantization of the progenitor halo masses. The algorithm 
also enables the merger process to be followed with high 
time resolution, as timesteps are not imposed on the tree 
but rather are controlled directly by the frequency of merg- 
ers. It is sim ilar in spirit to the method used by Kauffmann 
et al. ( |1993| ), but has several advantages, including smaller 
timesteps and not having to store large tables of progenitor 
distributions. Somerville & Kolatt (199i:) investigated a sim- 
ilar algorithm also based on equation (3.1), which they re- 
ferred to as binary mergers without accretion. They rejected 
that algorithm as it over-predicted the number of massive 
halos at high redshift, and instead opted for a more elaborate 
algorithm which t hey c ompared with N-body simulations in 
Somerville et al. (EOOOl). Our algorithm, which we first used 
in Baugh et al. ( |199^ ), differs from the one they rejected 
in two important respects. First, we explicitly account for 
accretion of o bjects below the mass resolution using the ex- 
pression (3.5). Second, we make the rather subtle choice of 
selecting the first progenitor mass. Mi, only in t he range 
Mrea < Ml < M/2 of the distribution defined by We 
have found, that together, these choices produce an algo- 
rithm that successfully produces distributions of progenitors 
which, on average, agree quite accurately with the analytic 
expressions given by the extended Press-Schechter theory. 
Moreover, statistics that are not predicted by the extended 
Press-Schechter theory, such as the frequency distribution of 
progenitors of a given mass (the extended Press-Schechter 
theory only predicts the mean of the distribution), we find to 
be in excellent agreement with the same statistics extracted 
from N-body simulations. This detailed investigation of the 
behaviour of the algorithm will be presented in a future pa- 
per (Lacey & Cole in preparation). 



3.1.2 Utilizing the Merger Trees 

Although the merger trees described above have very high 
time resolution, the nature of the galaxy formation rules that 
we implement below require placing the merger tree onto a 
predefined grid of timesteps. The original binary merger tree 
is used to find which halos exist at each timestep and to iden- 
tify which of them merge together between timesteps. As a 
consequence, mergers that are actually rapid, consecutive 
binary mergers in the original tree will appear as simulta- 
neous, multiple mergers in the discretized tree. The loss of 
information involved is not significant since, in reality, merg- 
ers are not instantaneous events and our discrete timesteps 
are typically much smaller than the dynamical timescales 
of the merging halos. Each merger tree thus starts from a 
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single halo of a specified mass M at z = 2haio, where Zhaio is 
the redshift for which we want to calculate the galaxy prop- 
erties. It extends up to some earlier redshift Zstart > Zhaio, 
where the tree has split into many branches. The generation 
of the halo merger tree proceeds backwards in time, starting 
from the trunk at z = Zhaio, but the calculation of galaxy for- 
mation and evolution through successive halo mergers pro- 
ceeds forwards in time, moving down from the top of the 
tree. The appropriate grid of A'stcps timesteps, the starting 
redshift Zstart and also the mass resolution, Mrcs, depend 
on the problem of interest. The timesteps can, for example, 
be chosen to be uniform either in time or in the logarithm 
of the cosmological expansion factor. For the models pre- 
sented here, we typically set Mrcs = 5 x Mq and use 
100 timesteps logarithmically spaced in expansion factor be- 
tween 2 = and 7. In our models, stellar feedback prevents 
significant amounts of star formation occurring in halos of 
very low circular velocity and so provided M^as is sufficiently 
low, the model results are not afi'ected by its value. 

The second way in which we manipulate the merger 
trees before applying our galaxy formation rules is by chop- 
ping each tree into branches that define the formation time 
and lifetime of each halo. So far we have done this using 
a simple algorithm. We start, at the first timestep, at the 
top of the tree (corresponding to the lowest mass halos in 
the merging hierarchy), and define each halo present as a 
new halo that formed at that timestep. We then follow each 
of these halos through their subsequent mergers until they 
have become part of a halo with mass greater than /form 
times the original mass. We normally set /form = 2. This 
point defines the end of the original halo's lifetime. Con- 
sistent with this definition, the point at which a new halo 
life begins is defined by the point when mergers produce 
a halo whose mass exceeds /form times the formation mass 
of its largest progenitor. When applying the galaxy forma- 
tion rules detailed in the following sections, we always treat 
the halos as if they retained, throughout their lifetime, the 
mass and other properties, (mean density, angular momen- 
tum, etc.) with which they formed. The mass accreted prior 
to the final merger, which can, in extreme cases, be as large 
as /form — 1 times the original halo mass, is effectively treated 
as if it were all accreted at the end of the halo's lifetime. We 
have chosen /form = 2 for consistency with our earlier work 



(Cole et al. 1994) in which a factor of two was built into the 
"block model" that we used to generate merger trees. With 
this choice the two models produce near identical results 
when the content and parameters of the galaxy formation 
model are the same. In Table 3 of Section 7.1, we include 
one variant model with /form = 1.5 which demonstrates that 
the model is not very sensitive to the choice of this parame- 
ter. This is natural as most halos end their lives when they 
are accreted onto much more massive halos and thus their 
lifetimes are robust to the choice of /form. 

In order to investigate the statistical properties of 
galaxy populations, we generate a set of merger trees start- 
ing from a grid of parent halo masses specified at some red- 
shift, Zhafo- For each given halo mass, we generate many 
realizations of the merger tree. The resulting model galaxies 
can then be sampled, taking account of the abundance of 
the parent halos at z = Zhafo, to construct galaxy catalogues 
with any desired selection criteria such as an absolute or 
apparent magnitude limit. Alternatively, properties such as 



the galaxy luminosity function or number counts can be es- 
timated directly by a weighted sum over the model galaxies. 
We have used the Press-Schechter mass function to estimate 
the halo abundance, but it is well known that this formula 
overestimates th e abu ndance of M-^ objec ts somewhat (e.g. 
Efstathiou et al. |198^ ; Lacey & Cole |1994| ). Rece ntly, Sheth, 
Mo & Tormen (200C) and Jenkins et al. (pOOCl) have pre- 



sented fitting formulae that match the results of N-body 
simulations to high accuracy. In future it will be preferable 
to use these formulae, but here we simply note that adopt- 
ing the Jenkins et al. (2000) mass function would make little 
difference to our model predictions. 



3.2 Halo Properties 

In order to calculate the properties of the galaxies that form 
within the dark matter halos produced by the merger tree, 
we need a model for the internal structure of the halos. This 
must specify the halo rotation velocity required to calculate 
the angular momentum of the gas that cools to form disks, 
and the halo density profile required to calculate the sizes 
and rotation speeds of the galaxies. 

The properties of dark matter halos formed in cosmolog- 
ical, coUisionless, N-body simulations have been extensively 
1985| , 



studi ed (e.g. Fren k et a l 



1987 ; Warren et al. 



White 



1995a, 



1992 



1996 



1988 



Barri es & Efstathiou 
1991: ; Navarro, Frenk 



1997 



Cole & Lacey . 

Moore et al. [I999b| , Jing |2000| ) 



The models detailed below are designed to be consistent 
with the results from these simulations. 



3.2.1 Spin Distribution 

Dark matter halos gain angular momentum from tidal 
torques operating during their formation. The magnitude 
of this angular momentum is conventionally quantified by 
the dimensionless spin parameter 



At 



Jh\Eh\ 



1/2 



(3.6) 



where Mh, Jh and Eh are the total mass, angular momen- 
tum and energy of the halo . The distributions of Ah found 
various N-body studies (^ arnes & Efstathiou 19871: 



tathiou et al. 1985 



Warren et al. 1992 



Efs- 



|Colc fc Lac~199e 



Lemson & Kauffmann 1999) agree very well with one an- 



other. They depend only very weakly on halo mass and on 
the form of the initial spectrum of density fluc tuati ons. 

A good fit to the results of Cole & Lacey (199£) is pro- 
vided by the log-normal distribution. 



P(AH)dAH = 



1 



27r(JA 



exp 



(In A - In Amcd)" 
2al 



Ah ' 



(3.7) 



with Amcd = 0.039 and a\ = 0.53. This fit was obtained 
specifically for halos with M* < Mh < 2M* in the case 
of an n = —2 power spectrum, which is the most relevant 
for CDM models on galaxy scales, but we stress that the 
fit parameters depend only very weakly on mass and on 
the slope of the power spectrum. For example, this fit also 
reproduces quite accuratel y the distribution plotted in Fig. 4 
of Lemson & Kauffmann (1999), which is for galactic mass 
halos in a rCDM simulation. We use this distribution to 
assign, at random, a value of Ah to each newly formed halo. 
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Note that we do not take account of a possible correlation 
between the angular momenta of merging halos. It would 
be necessary to do this if one wanted to follow the angular 
momenta of galaxy merger products, but we currently do 
not attempt this. 



aNFW, reduced by a factor of 2/3. Such a change has only a 
relatively small effect on the galaxy properties that we exam- 
ine below. The largest changes are to the disk scalelengths, 
which decrease by 10%, and to the disk circular velocities, 
which increase by 7.5%. 



3.2.2 Halo density Profile 

Our standard choice is to mo del the dark matter density 
profile using the NFW model (Navarro et al. 1995a): 



p{r) 



1 



^virpcrit 

/(tSNPw) ''/'■vir ('■/r-vir + flNFw)^ 



[r < rvir), (3.8) 



with /(aNFw) = ln(l + I/onew) — 1/(1 + tSNPw), truncated 
at the virial radius, r^ir. We define the virial radius as the ra- 
dius at which the mean interior density equals Avir times the 
critical density, pcrit ~ 3H'^ /{SivG). Here the virial overden- 
sity, Avir, is defined by the spherical collapse model which 
yields Avir = 178 for f^o = 1. Expressions for Avir in low 
f2o, open and fiat, universes can be foun d in th e appendices 
of Lacey & Cole ( |l993| ) and Eke et al. ( [l996| ) respectively. 
Confirmation that this definition of the virial radius is phys - 
ically sensible is provided by Fig. 13 of Cole fc Lacey (1996) 
and Fig. 10 of Eke, Navarro & Frenk (|l99^). These show 



that on average the transition between dynamical equilib- 
rium and the surrounding infall occurs close to this radius. 
The NFW profile has one free parameter, onkw, which is a 
scalelength measured in units of the virial radius. Allowing 
this one parameter (equivalent to the inverse of the concen- 
tration in the terminology of Navarro et al. ) to vary, the 
density profile accurately fits the profiles of isolated halos 
grown in cosmological N-body simulations for a wide rang e 
of masses and initial conditions (Navarro et al. 1996 1997 ), 
including simulations th at contain adiabatic gas as well a s 
collisionless dark matter ( Eke et al. 199^ ; Frenk et al. 1999 ) . 
Furthermore, there is a correlation between the best fit value 
of flNFW and halo mass. This can be understood in terms of 
how the typical forrn ation time of a halo depends on mass 
( Cole fc Lacey 199t ) . A simple analytic model for this re- 
lation has been presented in the appendix of Navarro et al. 
(1997), and it is this th at we use to set the values of ankw for 
our halos. Jing (2000) and Bullock et al. (199£) have found 
that there is considerable scatter about the mean correla- 
tion, which is presumably related to the differing dynamical 
states and formation histories of the halos. We do not take 
this into account, but we note that simply including this 
scatter by randomly perturbing the ankw values has little 
effect of the resulting distributions of galaxy properties. For 
instance, the distribution of galaxy sizes is already broad as 
a result of its dependence on the very broad distribution of 
halo angular momenta. 



Subsequently to the Navarro et al. (1997) paper, there 
has been some debate as to the accuracy with which the 
NFW profile fits the very central regions of dark matt er ha- 
los simulated a t ver y high resolution (Moore et al. 19 
Kravtsov et al. 1998). When a concensus is reached it may 



be possible to improve this aspect of our modelling by adopt- 
ing a slightly modified density profi le. We note that the most 
recent simulations by Moore et al. ( 1999a ) yield density pro- 
files which are sl ightly more centrally concentrated than the 
Navarro et al. (1997) result. To an accuracy of 20% they 
can be fitted by NFW profiles, but with the scalelengths. 



3.2.3 Halo Rotation Velocity 

To compute the angular momentum of that fraction of the 
halo gas that cools and is involved in forming a galaxy, we 
need a model of the rotational structure of the halo. We as- 
sume that the mean rotational velocity, Kot, of concentric 
shells of material is constant with radius and always aligned 
in the same direction. This simple description is broadly con- 
sistent with the behaviour seen in the simulations of War- 



ren et al. (1992) and Cole fc Lacey (1996). The appropriate 
value of Viot can be related to the halo spin parameter. Ah, 
by evaluating, for the adopted ha lo m odel, the quantities 
defining Jh and Eh in equation (3.6). This calculation is 
described in Appendix ^ We obtain 



Vrot ~ A(aNFw)AHVH, 



(3.9) 



where Vn = (GM/rvir)^^^ is the circular velocity of the halo 
at the virial radius. The dimensionless coefficient ^(ankw) is 
a weak function of flNFW, varying from A ~ 3.9 for anew = 
0.01 to 4.5 for onfw = 0.3. 

Our code allows us to explore the effects of using alter- 
native dark matter density profiles. In particular, we have 
included the case of a singular isothermal density profile, 
p{r) oc r~^, and a non-singular isothermal density profile, 
p{r) oc l/((r/rvir)^ + a^). We find A = 8\/2/7r ^ 3.6 for the 
singular isothermal sphere (see Appendix M). The value of A 
decreases very slowly as a core radius is introduced, falling 
to A^ 3.4 for a = 0.3. 

Our model of th e dist ribution of hot gas in the halo is 
described in Section |4.1.l| . As the hot gas is less centrally 
concentrated than the dark matter, if we were to assume 
they had identical rotation velocities this would result in the 
gas having a slightly higher mean specific angular momen- 
tum than the dark matter. We therefore take the rotation 
velocity of the gas to also be constant with radius, but with 
a value Vf^^ such the gas and dark matter have the same 
mean specific angular momentum within the virial radius. 
This simple model seems to be in reasonable accord with 
the properties of clusters in the high-resolution, gas-dynamic 
simulations of Eke et al. (1998, Eke private communication). 



4 FORMATION OF DISKS AND SPHEROIDS 

In this section we describe how disks and spheroids form, 
how we model star formation, feedback and chemical evolu- 
tion, and how we calculate galaxy sizes. 



4.1 Disk Formation 

We assume that disks form by cooling of gas initially in 
the halo. Tidal torques impart angular momentum to all 
material in the halo, including the gas, so that gas which 
has cooled and lost its pressure support will naturally settle 
into a disk. Below, we detail how we compute the mass of 
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the forming disk based on the radiative cooling rate of the 
halo gas, and how we compute its angular momentum. 



4.1.1 Hot Gas Distribution 

Diffuse gas which is not part of galaxies is assumed to be 
shock-heated during halo collapse and merging events. We 
will refer to this halo gas as "hot", to distinguish it from 
the gas in galaxies, which we call "cold". To calculate how 
much of this hot gas cools to form a disk, we need to know its 
initial temperature and density profile. In contrast to most 
previous work, we will not assume that the hot gas has the 
same density profile as the dark matter. 

High-resolution hydrodynamical si mulatio ns of the for- 
mation o f gala xy clusters ( Navarr o et al. 1995a ; Eke, Navarro 
& Frenk [L998| , Frenk et al. ^99§) show that, in the absence 
of radiative cooling, the resulting dark matter distribution is 
well-modelled by an NFW profile, but that the shock-heated 
gas is less centrally con centrated. The gas distribution i s 
well fit by the /3-model (Gavaliere & Fusco-Femiano 1976), 
pgas(r) oc (r^ + rcore)"^*"'''^, traditionally used to model 



the hot X-ray emi tting gas in galaxy clusters. The simula- 
tions of Eke et al. (1998), which span a narrow range of halo 
mass in an Qq ~ 0.3, Ao = 0.7 cosmology, indicate that 
the typical cluster gas profile is accurately described by a 
/3-model with « 2/3 and rcorc/j'NFW ~ 1/3. Here, rNPw 
is the NFW scalelength, equal to aNPWfvir, and so for these 
clusters rcore/j'vir ~ 1/20. A similar result was fo und fo r 
clusters in an f2o = 1 cosmology by Navarro et al. (1995a). 



In both cases, the simulations produce cluster gas temper- 
ature profiles that vary slowly with radius, consistent with 
hydrostatic equilibrium. The mean temperature of the gas 
is close to the virial temperature, defined by 



_ 1 tmvR 2 



(4.1) 



where mn is the mass of the hydrogen atom and ^ the mean 
molecular mass. 

Motivated by these simulation results, we assume that 
any diffuse gas present in the progenitors of a forming halo 
is shock-heated during the halo formation process and then 
settles into a spherical distribution with density profile, 



(r) oc l/(r -f Tcoro)- 



(4.2) 



For simplicity, we assume the gas temperature to be con- 
stant and equal to the virial temperature, Tvir. The effect 
this assumption has on the cooling radii and masses, com- 
puted below, is generally very small, as the cooling time of 
the gas depends more strongly on density than temperature 
and the density gradient is typically much larger than the 
temperature gradient. Guided also by the numerical sim- 
ulations, we assume that, for the first generation of halos, 
f'core = rNFw/3. Howcver, this result is for simulations which 
do not include radiative cooling, and we expect this rela- 
tionship to be modified for halos formed from progenitors in 
which gas has already been removed by cooling. The gas that 
is able to cool most efficiently in any halo is the densest gas 
with the lowest entropy. Thus, the remaining gas involved in 
the formation of a new halo will have a higher minimum en- 
tropy than if cool ing h ad not occured. The a nalytic work of 
Evrar d & Henry ( |l99l| ), Kay & Bower ( |l999| ) and Wu et al. 
(2000) suggests that increasing the minimum entropy of the 



halo gas has the effect of increasing its core radius. Further 
out, where cooling has had little effect, the gas properties 
will be less affected and, in particular, the pressure at the 
virial radius, which is ultimately maintained by shocks from 
infalling material, will remain unchanged. 

As a qualitative description of the behaviour described 
above we have constructed the following simple model. 
When a new halo is formed in a merger, if the hot gas frac- 
tion in the halo is less than the global value of f2b/f2o (in- 
dicating that some gas has already cooled), we increase the 
gas core radius, rcore, until we recover the same density at 
the virial radius that we would have obtained had no gas 
cooled. In principle, this ceases to be possible once the gas 
fraction is so low that even if it were placed in the halo 
at constant density, this density would be below the target 
value. To deal with this contingency, we set an upper limit of 
'"core = lOrvir, but iu practicc this extreme is rarely reached. 
The result of this procedure, when applied to the models 
discussed in Section |^, is that at high redshift the core radii 
start with values close to rcoro = rNPw/3, which for iso- 
lated bright galaxies, groups and clusters are approximately 
20,3Q,50/i^^ kpc respectively. As gas cools and galaxy for- 
mation proceeds, the core radii grow until, at the present 
day, the corresponding median core radii for newly formed 
halos are 85, 125 and 175/i~^ kpc. The distributions of core 
radii typically span a factor of two in scale. 

As alternatives to this standard description, our code 
also allows us to keep the core radius fixed, either as a fixed 
fraction of the virial radius or of the NFW scalelength, or 
even simply to assume that the gas traces the dark matter 
density profile. These options allow us to gauge directly the 
effects of our model assumptions. 

4.1.2 Cooling 

Assuming that the shock-heated halo gas is in coUisional 
ionization equilibrium, the cooling time, defined as the ratio 
of the thermal energy density to the cooling rate per unit 

volume, /9gasA(rgas, ^gas), is 



'^cool 



[r) 



2 firriH pgas(r)A(Tgas, Zg- 



(4.3) 



Here Pgas('") is the density of the gas at radius r, Tgas is 
the temperature and Zgas the metallicity. We use the cool- 
i ng fu nction A(rgas, ^gas) tabulated by Sutherland & Dopita 
( 1993| ). We estimate the amount of gas that has cooled by 



time t after the halo has formed by defining a cooling ra- 
dius, rcooi(i), at which TcooI = t. Note that for the purpose 
of computing this cooling radius, the gas density profile is 
kept fixed throughout the halo lifetime. 

The gas that cools is assumed to be accreted onto a 
disk at the centre of the halo. We estimate the time taken 
for this material to be accreted onto the disk as the free- 
fall time in the halo with the assumed density profile. Con- 
versely, we can define a free-fall radius rg{t) beyond which, 
at time t, material has not yet had sufficient time to fall 
into the central disk. Thus, to compute the mass that cools 
and is added to the disk in one timestep, At, we compute 
rmin(t) ~ min[rcooi, rff] at the begining and the end of the 
timestep, and set Mcooi At equal to the mass of hot gas orig- 
inally in the spherical shell defined by the two values of rmin. 
For one timestep, this defines the cooling rate, Mcooi, that 
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ente rs in to the differential equations (i.6) to (1.11) of Sec- 



tion 4.2 



describing the star formation, chemical enrichment 



and feedback. 



4-1.3 Angular momentum 

We assume that when the halo gas cools and collapses down 
to a disk, it conserves its angular momentum. Thus, the 
specific angular momentum of the material added to the 
disk by cooling since the formation of the halo is equal to 
that of the gas orig inally within rmin = min[rcooi, fff] ■ As 
described in Section 3.2.3 we take the rotation velocity of the 
hot halo gas, V^^^ , to be constant with radius, which implies 
that the specific angular momentum increases linearly with 
radius in the halo. 

The assumption of angular momentum conservation 
during the collapse is not a trivial one. In fact, numerical 
hydrodynamical simulations of galaxy formation including 
radiative cooling have, up to now, found that the cold gas 
loses most of its angular mome ntum (e.g. White & Navarro 
Navarro, Frenk & White |l995b[ Navarro & Steinmetz 
199£). However, these simulations have either not included 



star formation and feedback, or only included it in a very 
simple way which may not be accurate. In the absence of 
stellar feedback the gas distribution in a forming galactic 
halo is very clumpy. These clumps are efficient at losing 
angular momentum to the dark matter halo via dynami- 
cal fric tion. It is precisely this process that we model, in 
Section 



4.3.1 



to follow the merging of galaxies. However, if 
feedback keeps the gas that is not in galaxies diffuse, then 
the loss of angular momentum will be much reduce d. Thi s 
has been investigated by Weil, Eke, fc Efstathiou (1998), 
Sommer-La rsen, Gelato & Vedel ( ]199£| ) and Eke, Efstathiou 
& Wright (1999) who found that delaying the cooling of 
the gas considerably reduces the loss of angular momentum. 
We also note that strong angular momentum loss results in 
galaxy disk sizes much smaller than observed. In contrast, as 
we show later, our assumption of angular momentum con- 
servation leads to disk sizes very similar to observed values. 

In the following section we will introduce a model of 
stellar feedback whereby gas can be ejected from the disk. 
When this occurs, we assume that the specific angular mo- 
mentum of the remaining material is unaffected. In Section 
4.4 we give details of how we relate the size of the disk to 



its mass and angular momentum. 



4.2 Star Formation in Disks 

We now turn to the important process of star formation 
within disks. Star formation not only converts cold gas into 
luminous stars, but it also affects the physical state of the 
surrounding g SNe and young stars inject energy and 

metals back into the ISM. The energy that is released can 
be sufficient to drive gas and metals out of the galactic disk 
in the form of a hot wind. The removal of material from 
the disk acts as a feedback process which regulates the star 
formation rate. Also, the injected metals enrich both the 
cold star-forming gas and the surrounding diffuse hot halo 
gas. Enrich men t of the halo gas decreases the cooling time 
defined in (4.3), allowing more gas to cool at late times, 
while stellar enrichment affects the colour and luminosity of 



the stellar populations. Early semi-analytic models were un- 
able to include these effects accurately, as stellar population 
synthesis models with a wide range of metallicities were not 
available. Now that such models are widely available, simple 
chemical enrichment models have been included in several 



semi-analytic models, e.g. Kauffma nn ( 1996), Kauffmann & 
Chariot ( |l998a|' ), Guiderdoni et al. ( [1998D , and Somerville & 
Primack (11999). 



4-2.1 Chemical Enrichment and Feedback 

Our basic model of star formation assumes that stars are 
formed in the disk at a rate directly proportional to the 
mass of cold gas. Thus, the instantaneous star formation 
rate, i/>, is given by 



ii = Mcoid/r* 



(4.4) 



where the star formation timescale is r*. To model the feed- 
back effects of energy input from young stars and SNe into 
the gas, we assume that cold gas is reheated and ejected 
from the disk at a rate 



Me- 



(4.5) 



In general, both r^. and (3 are functions of the properties of 
the surrounding galaxy and halo. We will return later (Sec- 
tion 4.2.2) to the way in which we model these dependencies. 



The processes of gas cooling from the reservoir of hot 
halo gas and accreting onto the disk, star formation from 
the cold gas, and the reheating and ejection of gas all occur 
simultaneously. For each halo, we estimate the rate at which 
gas cools and is accreted by the central galaxyby comput- 
ing the cooling radius, as described in Section 4.1.2, at each 



discrete timestep at which the halo merger tree is stored. 
Within one of these discrete steps we approximate the cool- 
ing rate as a constant, Mcooi, and use a simple instantaneous 
recycling approximation t o model star formation, feedback 
and chemical enrichment (Tinsley 1980). Note that for satel- 
lite galaxies Mcooi = 0, as their hot gas is assumed to be 
stripped. Fig. ^ depicts the various channels by which mass 
and metals are transferred between the three phases. Note 
that we always compute Mcooi from the initial density pro- 
file of the hot gas and so we are implicitly assuming that gas 
reheated by SNe plays no role until it is incorporated into a 
new halo as a result of a merger. Under the instantaneous 
recycling approximation, the rate of fiow down each chan- 
nel is simply proportional to the instantaneous star forma- 
tion rate, ip, or the cooling rate, A'/coo1. The labels in Fig. ^ 
give the rates in terms of these quantities. The solid lines 
refer to total rates and the dashed lines to the metal com- 
ponent. Note that we have allowed for the possibility that 
some fraction of the metals produced by stars may be di- 
rectly transferred to the hot halo gas, but we have neglected 
the corresponding transfer of mass. This is a good approx- 
imation, since the directly ejected material would be very 
metal rich and the mass transferred by this route will al- 
ways be small compared to that transferred by reheating of 
the cold gas by SN feedback. 

In Fig. ^ and below, p denotes the yield (the fraction 
of mass converted into stars that is returned to the ISM in 
the form of metals), R the fraction of mass recycled by stars 
(winds and SNe), e the fraction of newly produced metals 
ejected directly from the stellar disk to the hot gas phase. 
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Figure 3. A schematic diagram showing the transfer of mass and metals between stars and the hot and cold gas phases during a single 
timestep. The solid lines indicate the routes and rates by which mass is transferred between the three reservoirs, while the dashed lines 
refer only to the exchange of metals. The instantaneous rate of star formation is ip and the cooling rate is Mcooi- The metallicities of 
the cold gas, stars and hot halo gas are ^coldi a^nd ^hot respectively. The yield of the assumed IMF is p and the parameters /3 and e 
describe the effect of SN feedback and the direct ejection of SN metals into the hot halo gas. 



.^coid the metallicity of the cold gas, and (3 the efhciency of 
stellar feedback. Each of the arrows in Fig. ^ gives rise to 
a term in the following differential equations that describe 
the evolution of the mass and metal content of the three 



reservons; 
M* 

Mhot = 

Afcold = 

+ 



(1 - R)^p 

-Mcool + Pip 

A'/cooi-(l-i? + /3)V 
(1 - R)Z,M^ 

- Mcool -^hot + (pe + (3Zcoid)tp 

Afcool.^hot 

- e) - (1 + /3 - i?)Zcoid)V', 



(4.6) 
(4.7) 
(4.8) 
(4.9) 
(4.10) 

(4.11) 



where Zcoid = A^coid Mucoid and Zhot = Mif„t/Mhot. The 
values of R and p in these equations are related to the IMF, 
as discussed in Section |3.2| . 

We assume that over one timestep the cooling rate, 
Mcool, and the metallicity of the hot gas, Zhot, can be taken 
to be constant. This set of first-order, coupled differential 
equations can be straightforwardly solved to give the change 
in mass and metal content of cold gas, hot gas and stars since 
the start of the timestep (Appendix ^) . The model is quite 
flexible: its behaviour is determined by specifying how the 
functions r*, f3 and e depend on the properties of the galaxy 
and its surrounding halo. We note that compared to the 
simple, "closed-box" chemical enrichment model, the yield 
is modified by the metal ejection and feedback to produce 
an effective yield Pcff = (1 ~ e)p/(l ^ R + (equation B9), 
which is therefore a function of the potential-well depth of 



the galaxy. The evolution of the stellar metallicity differs 
from the closed-box model because it is affected by both the 
ejection of reheated gas and the accretion of cold gas and 
associated metals. 



4-2.2 Star Formation Law and Feedback Parameterization 

In our previous work (e.g. [Cole et al. 1994 ) , we specified the 
star formation timescale and feedback efficiency in terms of 
the circular velocity of the halo in which each galaxy formed, 
Vh. T he relations we adopted were 



'(Vh/300 kms" 



and 



(4.12) 



(4.13) 



The parameter r°', we treated as a free parameter, while the 
other three parameters, ai, V]^^^ and adoti we constrained 
by comparing our models to the num erica l simulations of 
galaxy formation of Navarro & White (1993). These simula- 
tions had only one free parameter, the fraction of SN energy 
injected as kinetic energy into the interstellar medium. In 
order to suppress the formation of low luminosity galaxies, 
and thus produce a galaxy luminosity function with a rea- 
sonably shallow faint end slope, as observed, we adopted 
a fiducial model with very strong feedback for low circular 
velocity halos, which we obtained by setting the parameter 
values = —1.5, Vi'^t = 140kms~^ and a'^^at = 5.5. 

The more detailed modelling that we now perform of the 
structure of our model galaxies allows us to specify the star 
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formation timescale and feedback efficiency more naturally 
in terms of the properties of the galaxy disk, namely its cir- 
cular velocity, Vdisk, and dynamical time rdisk = ?'disk/Vdisk- 
Vdisk and rdisk are both taken at the disk half-mass radius. 
The relations that we adopt are 



-r* = ^ Tdisk (Vdisk/200kms ^)"* 
and 

/?= (Klisk/Hot)-"""*, 



(4.14) 



(4.15) 



where e*, a* and ahot are dimensionless parameters, and the 
parameter, Vhot, has the dimensions of velocity. If a* = 0, 
then our star formation law, (4.14), simply gives a star for- 
mation timescale proportional to the galaxy dynamical time, 
broadly consistent with the observational data compiled by 
Kennicutt (199S). The inclusion of the velocity dependent 



term allows us to explore models that have a similar depen- 
dence on velocity as our p revio us, quite successful, model 
which had a^ = —1.5 in ( [l.l2| ). It should be noted that 
because the cold gas reservoir is depleted both by the for- 
mation of stars and by reheating due to SN feedback, the 
timescale on which the reservoir is depleted (in the ab- 
sence of any further gas cooling) is shorte r th an r ^. Iri Ap- 
pendix^ where the analytic solutions of ( [4.6| ) to (4.11) are 
discussed, it is shown that this timescale, which in turn de- 
termines the effective star formation timescale, is given by 
Toff = r*/(l — JJ-I-/3). The feedback equation is the same 
as we used previously, but now expressed in terms of the 
galaxy circular velocity rather than the halo circular veloc- 
ity. This is physically more realistic as it is the depth of the 
potential at the point where the stars are forming which is 
most relevant. To constrain these four parameters we now 
prefer to take a more empirical approach and use a wider 
range of observational data, rather than to fix the parame- 
ters to emulate one particular set of numerical simulations 
of galaxy formation, as we did before. 



4.3 Spheroid Formation 

In our model, the primary route by which bright elliptical 
galaxies and the bulge components of spiral galaxies form 
is through galaxy mergers. When dark matter halos merge, 
we assume that the most massive galaxy automatically be- 
comes the central galaxy in the new halo, while all the other 
galaxies become satellite galaxies orbiting within the dark 
matter halo. The orbits of these satellite galaxies will grad- 
ually decay as energy and angular momentum are lost via 
dynamical friction to the halo material. Thus, eventually, the 
satellite galaxies spiral in and merge with the central galaxy. 
We now describe how we estimate the times at which such 
galaxy-galaxy mergers occur and what they produce. 



4-3.1 Dynamical Friction 

When a new halo forms, we assume that each satellite galaxy 
enters the halo on a random orbit. The most massive pre- 
existing galaxy on the other hand is assumed to become 
the central galaxy in the new halo, where it will act as the 
focus for any gas that may cool within the new halo. The 
time for a satellite's orbit to decay due to the effects of 
dynamical friction depends on the initial energy and angular 



momentum of the orbit. Lacey & Cole ( 1993| ) estimated the 
time for an orbit to decay in an isothermal halo, based on the 
standard Chandrasekhar formula for the dynamical friction. 
Their formula (B4) can be written in the form, 



'^mrg — fdf ©orbit "^dyn 



0.3722 



ln(Ac 



(4.16) 



Here, Mh is the mass of the halo in which the satellite orbits, 
and we take Msat to be the mass of the satellite galaxy 
includin g the mass of the dark matter halo in which it 



formed (Navarro et al. 1995c). Note that we deliberately 
count the mass of the satellite's halo in the definition of 
both Maat and Mh. The Coulomb logarithm, we take to be 
In(Acouiomb) ~ \n{M-ii/Mazx). The dynamical time of the 
new halo is rdyn = Trrvir/Vk, defined equivalently as either 
the half period of a circular orbit at the virial radius, or as 
(Gpvir)"^'^^, where pvir is the mean density within the virial 
radius, or, for an isothermal sphere, as the full orbital period 
of a circular orbit at the half-mass radius. 

The dependence of this merger timescale, Tmrg, on the 
orbital parameters is contained in the factor Oorbit, defined 
as, 



(4.17) 



where E and J are the initial energy and angular momentum 
of the satellite's orbit, and rc{E) and Jc{E) are the radius 
and angular momentum of a circular orbit with the same 
energy as that of the satellite. The power-law dependence 
on the circularity, J/Jc{E), is an accurate fit to the result 
of numerical integration of the orbit-averaged equations de- 
scribing the ef fect of dynamical fi ction in the range 0.01 < 
J/Jc{E) < 1 dLacey fc Cole 1993| ). The distribution of ini- 
tial orbital parameters of infalling satellites in cosmological 
N-body simulations has been studied by Tormen (1997). We 
find from his results that the distribution of Oorbit is well 
modelled by a log normal with (logj^Q Oorbit) ~ —0.14 and 



2x1/2 



0.26. 



dispersion ((log^o Oorbit - (logio Oorbit)) ) 

The merger timescale computed in this manner is based 
on several approximations, e.g. treating the satellite as a 
point mass with mass equal to the sum of the galaxy mass 
plus that of its original dark matter halo. We therefore allow 
ourselves some freedom by inserting the dimensionless pa- 
rameter /df , which is greater than unity if the infalling satel- 
lite's halo is efficiently stripped off early on. We note that 
recent analyt ical a nd numerical investigations by va n den 
Bosch et al. ( |l999| ) and Colpi, Mayer & Governato ( |199E| ) 
suggest a weaker dependence of the merger timescale o n the 
orbital circularity, with the exponent 0.78 in equation (4.17) 



being replaced by a value of 0.4 or 0.5, but these results were 
also derived using a somewhat different halo density profile 
from the s ingular isothermal sphere assumed by Lacey & 
Cole (1993). I n this work we h ave retained the model defined 
by equations (4.16) and (4.17), but we note that it may soon 
be possible to have a fully specified and calibrated model for 
dynamical friction-driven mergers. 

The procedure for determining the fate of satellite 
galaxies within dark matter halos is straightforward. When a 
new halo forms, each of the satellite galaxies that it contains 
is assigned a random value of Oorbit according to the log- 
normal distribution described above. Then, for each satel- 
lite, we compute Tmrg from equation (4.16). The satellite is 



assumed to merge with the central galaxy after this time in- 
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terval has elapsed, provided this occurs during the lifetime 
of the halo, i.e. before the halo has merged to become part 
of a much larger system. Satellites that do not merge are as- 
signed a new random value of Oorbit when the halo in which 
they reside is incorporated into a new, more massive halo. 



4-3.2 Galaxy Mergers and Bursts 

Our method for modelling galaxy mergers produces, at each 
timestep, a list of satellite galaxies which merge with the 
central galaxy in each halo. If our grid of timesteps were 
sufficiently fine, then these lists would always contain just 
one or zero satellite galaxies, but in practice there is often 
one large satellite and several smaller satellites merging with 
the central galaxy at a single timestep. We deal with this by 
ranking the merging satellites by mass and then, starting 
with the most massive one, merge them sequentially with 
the central galaxy. 

The outcome of each merger depends on the ratio of 
the mass of the merging satellite, Maat , to that of the cen- 
tral galaxy, Mccn, a nd ha s been studie d rece ntly by Walker, 
Mihos & Hernquist (1996) and Barnes (1998), using numeri- 
cal simulations. As a simplified description of the outcome of 
these merge rs, we adopt the pr escripti on used in Kauffmann 
et al. (|l99j ) and Baugh et al. ( |l996a|) : 



in the context of the hiera rchical formation of galaxies by 
Mo, Mao & White (|l998a|). Our disk model is similar to 



a) If the mass ratio of merging galaxies, defined in terms 
of stars and cold gas only, is Msax/Mcan > .foUip, then the 
merger is said to be "violent" or "major", and a single bulge 
or elliptical galaxy is produced. Any gas present in the disks 
of the merging galaxies is converted into stars in a burst. 
We use the standard star formation and feedback rules, but 
now based on the circular velocity and dynamical time of the 
spheroid that is formed rather than the disk, and with a very 
much shorter timescale, similar to the dynamical timescale 
of the spheroid. 

b) Alternatively, if Afaat/Mcen < /eiiip, then the merger is 
classed as "minor", and, unless explicity stated otherwise, 
the stars of the accreted satellite are added to the bulge of 
the central galaxy, while any accreted gas is added to the 
main gas disk without changing the disk's specific angular 
momentum. 

The merger simulations mentioned above have not been run 
for a wide enough range of initial conditions to determine 
/ciiip exactly, but suggest a value in the range 0.3 /oiiip ^ 
1. The way in which we calculate the size of the spher oid 
which forms from a merger is described in Section 4.4.2, In 
the case of minor mergers, we also have the option of adding 
the accreted stars to the disk of the central galaxy. If we do 
this, we assume that the specific angular momentum of the 
disk is unchanged by the accretion. 



4-3.3 Disk Instabilities 

An issue we have not yet addressed is whether the disks 
in our model galaxies are dynamically stable. In particular, 
strongly self-gravitating disks are likely to be unstab le to the 
formation of a bar (e.g. Bin ney fc Tremaine 



tathiou. L ake fc Negroponte 
& Tohhne 



1995; Sellwood 



1987, ^6; Efs- 



1982 ; Christodoulou, Shlosman 
1999|; Syer, Mao fc Mo |l999| ). Re- 



theirs, except that we explicitly follow the formation and 
structure of a bulge component and, more importantly, we 
follow the complete merging history of both the bulge and 
the dis k. The stability criterion considered by Mo, Mao fc 
White (1998a) is based on the quantity: 



(4.18) 



(GMdisk/rdisk)i/2- 

According to Efstathiou, Lake & Negroponte (1982), for 
disks to be stable requires > 1.1. In the original for- 
mulation, Knax was the rotation velocity at the maximum 
of the rotation curve, but in our models we use instead the 
circular velocity at the disk half-mass radius. 

We have an option in our code to include the effect 
of such disk instabilities on g alaxy evolution. In that case, 
we check the criterion (4.18) for each galaxy disk at each 



timestep. If at any point a disk is unstable according to 
this condition, we assume that the instability results in the 
stellar disk evolving into a bar and then into a spheroid 



( [Combes et al. 1990| ; |Combcs 1999| ). We also assume that 



bar instability causes any gas present in the disk to undergo 
a burst of star formation subject to our standard feedback 
prescription. 

We do not include the effects of disk instabilities in our 
reference model. We brieffy present the effect it has on the 
distribution of disk sc alele ngt hs a nd the morphological mix 
and 



of galaxies in Sections 7.3 



7.4 



cently, the incidence of unstable disks has been considered 



but we postpone to a fu- 
ture paper a more detailed exploration of their consequences. 



4.4 Galaxy Sizes 

The two basic principles upon which we base our estimates 
of galaxy sizes are: 

i) the size of a disk is determined by centrifugal equilib- 
rium and conservation of angular momentum 

and 

ii) the size of a stellar spheroidal remnant produced by 
mergers or disk instability is determined by virial equilib- 
rium and energy conservation. 

The application of these simple principles is complicated by 
the gravitational interaction of the galaxy disk, spheroid and 
surrounding dark matter halo. Because of this, to determine 
either the disk or bulge radius, we must solve for the simul- 
taneous dynamical equilibrium of all three components. We 
use the following approach: 

a) The disk is assumed to have an exponential surface den- 
sity profile, with half-mass radius rdisk- 

b) The spheroid is assumed to follow an r^^'^ law in pro- 
jection, with half-mass radius (in 3D) rbuigc. 

c) The dark halo has a specified initial density profile 
(NFW in the standard case), but this is spherically deformed 
in response to the gravity of the disk and spheroid. 

d) The mass distribution in the halo and the lengthscales 
of the disk and bulge are assumed to adjust adiabatically 
in response to each other: for the disk, we assume that the 
total angular momentum is conserved; for the halo, we as- 
sume that rVc{r) is conserved for each spherical shell; for 
the spheroid, we assume that rV'c(r) is conserved at rbuige- 
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The task is then to solve for rdisk, ^buige and the deformed 
halo profile in dynamical equilibrium, subject to these con- 
straints. The method is described in detail in Appendix ^ 
This adiabatic invariance method for calculating the re- 
sponse of a halo or spheroid to the dis k was originally devel- 
ope d and applied by Barnes & White (1984), Blumenthal et 
al. (Il986|) and Ryden & Gunn (|1987|). 



4.4-1 Disk Sizes 

As already stated, the size of a disk is basically deter- 
mined by the angular momentum of the halo gas which 
cools to form it. Many previous papers have used a ver- 
sion of the following argument: if the dark halo and the 
gas it contains are modelled as a singular isot herm al sphere 
(p oc r~^, then, from the results of Section 3.2.3 and Ap- 
pendix Wi the mean specific angular momentum of the gas 



which cools is 



o T^COol 



= y2AHrcooiVH. On the 



other hand, if the self-gravity of the disk is also neglected, 
it rotates at constant circular velocity Vk, and so has mean 
specific angular momentum jdisk = '2huVu, for an expo- 
nential disk with scalelength ho- Equating jcooi and jdisk 
gives rdisk ~ l-68hn = I-WXhTj^i- This simple relation 
was originally derived by Fall ( [l983| ). It was used to calcu- 
late disk sizes i n gala xy formation models (wit h a fi xed Ah) 
by Lac ey et al. ( |l993| ), Kauff'mann & Cha riot ([L994| ), Kauff- 
mann (1996) and Somerville & Primack (1999). In this pa- 
per, we improve on this simple calculation by including (a) 
non-isothermal halo profiles for the dark matter and gas; 
(b) an initial distribution of Ah; (c) disk self-gravity and (d) 
gravity of the halo and spheroid, and their contraction in 
response to the disk. Most of these improvements were also 
i nclude d in the work on disk sizes by Mo, Mao & White 
(1998a), using similar techniques to those used here. How- 
ever, their work did not include a physical model for galaxy 
formation, so that they were forced to treat the disk-to-halo 
mass ratio, the disk-to-halo angular momentum ratio, and 
the disk M/L ratio as free parameters. If for a given halo we 
adopt the same disk angular momentum and mass, then our 
model produces disk scalesizes that agree very accurately 
with the Mo, Mao & White model. 



4-4-2 Sizes of Spheroids Formed by Mergers 

Spheroids can form either in major mergers (when any pre- 
existing disks are destroyed) or in minor mergers (when the 
disk of the larger galaxy survives). To estimate the size of 
the spheroid formed, we assume that the merging compo- 
nents spiral together under the action of dynamical friction 
until their separation equals the sum of their half-mass radii. 
At this point, we assume that the systems merge together, 
and we use energy conservation and the virial theorem to 
compute the size of the remnant. These considerations lead 
to: 



Mt 

ri r2 



'I _^ /orbit Ml Ma 



ri -I- r2 



(4.19) 



to 



(Ml + Ma) 

^new 

which relates the half-mass radius of the remnant, 
the masses. Mi and Ma, and half-mass radii, ri and ra, of 
the merging components. Defining Mi > Ma, Mi is the total 
galaxy mass for a major merger and the bulge mass for a 
minor merger, while Ma is the total galaxy mass for a major 



merger and the total stellar mass of galaxy 2 for a minor 
merger. The masses Mi and Ma include contributions from 
the respective dark matter halos, which are taken to be twice 
the halo mass within the half-mass radii ri or ra. 

The form factor, c, and the constant, /orbit, are related 
to the gravitational self-binding energy of each galaxy. 



-Ebin 



and their mutual orbital energy. 



/orbit GMiMa 



(4.20) 



(4.21) 



2 ri+r2 

at the point at which the merger occurs. The value of c 
depends weakly on the density profile of the galaxy; c = 0.49 
for an exponential disk and c = 0.45 for an r^-'^-law spheroid. 
For simplicity, we adopt c — 0.5. For the orbital energy, we 
adopt /orbit = 1-0, which corresponds to the orbital energy 
of two point masses in a circular orbit with separation ri -|- 
ra. These assumptions lead to the result that, for a merger 
of two identical, equal mass galaxies, the half-mass radius 
of the remnant increases by a factor rncw/ri — 4/3, which 
agrees reasonably well with the fact or of 1.42 found in the 
simulated galaxy mergers of Barnes (1992). 

□w , we then 



Having solved equation (4.19) for rbu 



ulgG ■ 



adiabatically adjust the spheroid, disk (if any) and halo to 
find the new dynamical equilibrium, as described in Ap- 
pendix ^ Typically, this leads to little change in rbuigc, 
showing that our treatment of the dark matter during the 
merger is approximately self-consistent. 



4-4-3 Sizes of Spheroids formed by Disk Instabilities 

As mentioned in the previous section, our code has an option 
to form spheroids through bar instabilities in disks. In this 
case, we compute the size of the resulting spheroid using 
virial equilibrium and energy conservation in much the same 
way as for the spheroids produced by mergers. If the mass 
of the unstable disk is Mdisk, the mass of any pre-existing 
central stellar bulge is Mbuigo, and their respective half-mass 
radii are rdisk and rbuige, then we calculate the final bulge 
half- mass radius, rncw, from the relation 



CB (Mdisk -I- Mbulgc)^ 



CBMb.ige ^ cuMl.i, 



+ /i 



^bulgc ^disk 

-A-^bulgc-^disk 
lit — . 

^bulgc H~ '^disk 



(4.22) 



Here we adopt cd = 0.49 and cb = 0.45, the form factors 
appropriate for a n exp onential disk and r^/*-law spheroid 
respectively (see (4.20)). The last term represents the grav- 



itational interaction energy of the disk and bulge which is 
reasonably well approximated for a range of rbuigc /rdisk by 
this form with /int = 2.0. After we h ave c alculated the new 
spheroid radius rnew from equation ( 4.22 ), we adiabatically 



adjust the spheroid and halo to a new dynamical equilib- 
rium, as for the case of a spheroid formed by a merger. 



5 GALAXY LUMINOSITIES AND SPECTRA 

The aspects of the model described so far enable us to fol- 
low the star formation history, chemical enrichment and size 
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Figure 4. The evolution of the B-band luminosity and the B-V 
and V-K colours for a single age stellar population. The solid 
lines show results for a stellar population with a Salpeter IMF for 
three different metallicities. The middle curves are for solar metal- 
licity, Z = 0.02, and the lower and upper curves for Z = 0.008 
and 0.05 respectively. The absolute magnitudes are normalized to 
IMq of stars. The corresponding dashed curves show results as- 
suming the Kennicutt IMF. In this case, the luminosities in each 
band have been reduced by a factor of 1.69 to make the solar 
metallicity curves for the two IMFs cross at an age oi t = 15Gyr. 



evolution of each galaxy. In order to convert this information 
into observable properties, we must model the spectropho- 
tometric properties of the stars that are formed, and the 
effects of dust and ionized gas within each galaxy on the 
emerging integrated galaxy spectrum. The models we adopt 
for each of these processes are outlined below. 



5.1 Stellar Population Synthesis 

The tech nique of st ellar population synthesis, pioneered by 
Tinsley (1972; 1980 ) and developed by Guid erdon i & Rocca- 
Volmerange (1987 ), Br uzual & Chariot (1993), Bressan, 



Chiosi & Fagotto ( 1994 ) and others, enables the observable 



properties of a stellar population to be computed given an 
assumption about the stellar Initial Mass Function (IMF) 
and the star formation history. Th e latest models of Bruzual 
& Chariot (2000, in preparation) provide the spectral en- 
ergy distribution (SED), l\{t,Z), of a single population of 
stars formed at the same time with the same metallicity, 
as a function of both age, t, and metallicity, Z. These can 
be convolved with the star formation history of a galaxy to 
yield its SED: 



Lx{t) 



lx{t-t',Z{t')) V(t') dt', 



(5.1) 



where Z{t') is the metallicity of the stars forming at time t' , 
and 'ip{t') is the star formation rate at that time. In the case 
of a galaxy which formed by merging, we also sum the con- 
tributions to L\ from the different progenitor galaxies, each 
with their own star formation and chemical enrichment his- 
tory. In performing the convolution integral, we interpolate 
the grid of SEDs, lx{t, Z), provided by Bruzual & Chariot, to 
intermediate ages and metallicities using linear interpolation 
in t and log Z. Broad-band colours can then be extracted by 
integrating over these spectra weighted by the appropriate 
filter response function. 

In our models, we always assume that the IMF is uni- 
versal in time and space. Observationally, the IMF is best 
constrained in the solar neighbourhood. However, even here 
there is significant uncertainty arising mainly from ambi- 
guity in the past star formation history. Because of this, 
we consider two possible choices of IMF, the form pro- 
posed by Sa lpeter (1955) and the form proposed by Ken- 
nicutt (1983), both of which produce reasonable agreement 



with the solar neighbourhood data. The Salpeter IMF has 
dN/dlnm oc m'" with x = 1.35, while the Kennicutt IMF 
has X = 0.4 for m < Mq and x — 1.5 for m > Mq. In both 
cases, visible stars have O.IM© < m < 125M0. The Salpeter 
IMF has been widely used in modelling galaxy evolution be- 
cause of its simplicity and the fact that it fits the observa- 
tional data on high mass stars fairly well. However, there 
is now cons i derab le observational evidence, as reviewed by 
Scalo (1986, 199^ ), that the IMF slope at low masses is fiat- 



ter th a.n the Salpeter form. The "best" IMF propo sed by 
Scalo (1998), which supersedes that of Scalo (1986), is ac- 
tually quite close to that of Kennicutt ( |l98j ) . We therefore 
adopt the Kennicutt IMF as our standard choice. 

We also include in our assumed IMF brown dwarfs (m < 
O.IMq), which contribute mass but no light to the stellar 
population. The fraction of brown dwarfs is specified by the 
parameter T, defined as 



T = 



(mass in visible stars -I- brown dwarfs) 
(mass in visible stars) 



(5.2) 



at the time a stellar population forms, i.e. before taking 
account of the fraction R of the mass that is returned to 
the ISM by recycling. Thus, by definition, T > 1. The effect 
of including brown dwarfs is to reduce all stellar population 
luminosities by a factor 1/T. We will see in Section [TJ that 
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observational estimates of the mass-to-light ratios of stellar 
populations constrain viable models to have modest values 
of T in the range 1 < T ^ 2. 

The way in which the predicted luminosity and colour of 
a stellar population depend on age, metallicity and choice of 
IMF is illustrated in Fig. ^. A number of properties which af- 
fect the behaviour of our galaxy formation models are worth 
noting. The overall stellar mass-to-light ratio depends sig- 
nificantly on the choice of IMF. This dependence has been 
explicitly scaled out of the curves shown in Fig. ^ by re- 
ducing all the luminosities in the Kennicutt IMF case by a 
factor of T = 1.69, so as to force the solar metallicity curves 
for the two IMFs to agree at i = 15Gyr. The slope of the 
absolute magnitude versus time curve has some dependence 
on the choice of IMF. For example, as the age is reduced 
from 15Gyr to around 3Gyr, the stellar population with the 
Kennicutt IMF brightens more rapidly than that with the 
Salpeter IMF. The differe nce is even larger for an IMF such 
as the Miller-Scalo IMF ( [Miller fc Scalo 1979[ ) which con- 
tains a greater fraction of stars of a few solar masses. In 
spite of the dependence of luminosity on the IMF, the B-V 
and V-K colours, both as a function of age and metallicity, 
are quite insensitive to the choice of IMF. The colours do 
depend strongly on metallicity, with increasing amounts of 
metals producing redder stellar populations. A young stel- 
lar population is very blue but rapidly reddens during its 
first 5Gyr. At later times the dependence of colour on age 
is much weaker. 

There are many approximations and assumptions in- 
volved in constructing stellar population synthesis models 
such as those of Bruzual & Chariot. Because of this, the ac- 
curacy of the model predictions is difficult to quantify. This 
i ssue has been addressed by Chariot, Worthey & Bressan 
( 1996 ) by comparing model predictions from different codes 
and for varying sets of assumptions. Their results indicate 
that for the same choices of IMF and star formation history, 
the resulting broad-band colours can differ by a few tenths of 
a magnitude, and this could give rise to 20-30% uncertain- 
ties in either the inferred age or metallicity. While efforts 
have been made, and continue to be made, to improve these 
models, it should be noted that the uncertainties in the pop- 
ulation synthesis model are sufficiently small that, for our 
purposes, the dominant source of uncertainty in modelling 
galaxy formation is, instead, the choice of IMF and its as- 
sociated yield. 



5.2 Yield and Recycled Fraction 

There are two further quantities related to the IMF that 
significantly affect galaxy formation. These are the recy- 
cled fraction, R, an d the yield, p. They appeared in equa- 
tions (4.6) to (4.11) for the evolution of gas and star masses 
and metallicities. The material which goes to form massive 
stars is mostly released back into the ISM via stellar winds 
and SN explosions. The returned gas is an important source 
of fuel for forming further generations of stars. SN explo- 
sions also enrich the ISM with metals giving rise to sub- 
sequent generations of redder, more metal rich stars. The 
recycled fraction and the yield are defined so that for each 
mass, AM, formed in new stars (including brown dwarfs), 
a mass RAM is returned to the ISM, and a mass pAM of 
newly synthesized metals is released. These quantities are 



given respectively by integrating the total ejected mass and 
the ejected mass in newly synthesized metals over the IMF. 
We recall that in a closed-box model of chemical evolution, 
the mean metallicity of the stars asymptotes t o a v alue of 
p/{l ~ R) as the gas is exhausted (e.g. Tinsley 198C). 

The values of R and p for any specific IMF can be 
estimated from stellar evolution theory and models of su- 
pernova explosions. We have used two different compila- 
tions of stellar evolution calcu lations to set these parame- 
ters: (i) Renzini & Voli's (1981) for intermed iate m ass stars 
{1 < m <, 8Mq), and Woosley & Weaver's ( |l995| ) for mas- 
sive stars (m ^ 8Mq) which produ ce Ty pe II supernovae; 
and (ii) results from Marigo et al. (1l99q) for intermediate 
mass stars and from Portinari et al. ( [1998|) for massive stars. 
The more recent calculations in (ii) include the effects of con- 
vective overshooting and quiescent mass loss. However, they 
rely on the supernova calculations of Woosley & Weaver. 
The contribution of SNII to the yield is sensitive to the as- 
sumed explosion energy (Woosley & Weaver's cases A,B,C); 
we give below the corresponding range in p for case (i), 
but Portinari et al. only calculated results for Woosley & 
Weaver's case A (case C would give larger yields). Type I 
supernovae make only a small contribution to the net pro- 
duction of heavy elements, and are not included here. The 
results for solar metallicity are as follows: for the Kennicutt 
IMF, case (i) gives Ri = 0.42, pi = 0.013 - 0.023, and case 
(ii) gives Ri = 0.44, pi = 0.022; for the Salpeter IMF, case 
(i) gives Ri = 0.28, pi = 0.010 - 0.020, and case (ii) gives 
Ri = 0.30, pi = 0.018. These values assume that T = 1. 
If T > 1, the appropriate values become p — pi/T and 
R = Ri/T. As may be seen from these values, for a given 
IMF, the recycled fraction, R, is fairly accurately known, 
but the theoretically predicted yield, p, is uncertain by at 
least a factor of 2. In our modelling, we have chosen to set 
R according to the above estimates. However, as the yield 
is more uncertain, we use these estimates only as a guide 
for what is reasonable and instead rely on observed galaxy 
metallicities to constrain the value of p. 

5.3 Extinction by Dust 

Absorption of starlight by dust has a significant effect on 
the optical luminosities and colours of galaxies, and a large 
effect on the far-UV luminosities which are used as the main 
tracer of star formation rates at high redshift. We model the 
effects of dust in a phys ically self-consistent way, using the 
models of Ferrara et al. (1999). Ferrara et al. have calculated 
radiative transfer of starlight through dust, including both 
absorption and scattering by dust grains, for a realistic 3D 
distribution of stars and dust, giving the net attenuation of 
the galaxy luminosity as a function of wavelength and incli- 
nation angle. In their model, stars are distributed in both 
a bulge and a disk, and dust i s dist ributed smoothly in a 
disk. The bulge follows a Jaffe (1983) distribution (which is 
very similar to an r^'^'^-law) with projected half-light radius, 
Tc. The stars and dust in the disk both have radially and 
vertically exponential distributions. The dust is assumed to 
have the same radial scalelength, ha, as the stars, but its 
scaleheight, h^, is in general different. The total dust con- 
tent is parameterized by the central V-band optical depth, 
Tvo, defined as the extinction optical depth looking verti- 
cally through the whole disk at r = 0. The dust properties 
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are chosen to match observations of the extinction law and 
albedo of dust in either the Milky Way (MW) or Small Mag- 
ellanic Cloud (SMC). 

Ferrara et al. tabulate separately the attenuations of 
disk and bulge light, as functions of wavelength, A, in- 
clination, i, central optical depth, tvo, ratio of bulge-to- 
disk scalelengths, re/hn, and ratio of dust-to-stellar verti- 
cal scaleheights, /iz,dust/^z, stars- We choose a fixed value for 
'Jz, dust //iz, stars, and Calculate tvo and rc/hR for each galaxy 
directly from the output of our model. We assign our galax- 
ies random inclination angles, and then calculate the at- 
tenuation factors for the disk and bulge luminosities at the 
wavelengths of each of the filters (e.g. B, K) we are using 
by interpolating in the tables. 

We calculate tvo for our model galaxies by assuming 
that it scales as the dust mass per unit area which, in turn, 
is assumed to scale with the total mass of metals per uint 
area in the cold gas: 



Mdust A'fcold^cold 

TVO OC — CX 2 ■ 

''disk ''disk 



(5.3) 



The metallicity Zcoid is obtained from our chemical evolution 



calculation. We normalize equation (5.3) by assuming that 
gas with solar metallicity, Z — .02, h as the local ISM dust- 
to-gas ratio. Savage & Mathis ( 1979| ) find Av/Nh = 3.3 x 
10~'^^magcm^ for the local ratio of V-band extinction, Av, 
to hydrogen column density A'^e. This then implies 



TVO = 0.043 



Mcoid/(27r/i| 



M, 



opc- 



0( 



■^cold 

0.02 



(5.4) 



Our standard choice is to use an MW extinction curve 
and to assume ftz,dust/^z, stars = 1. We have investigated 
variations in /iz,dust/ftz, stars over the range 0.4 to 2.5, and 
find that most results are very insensitive to this. Most re- 
sults also do not change significantly if an SMC rather than 
an MW extinction curve is used. The SMC and MW extinc- 
tion curves differ significantly only in the far-UV, but even 
here, the effects on our results are fairly small, as the net 
attenuation of galaxy light calculated using these radiative 
transfer models has a much weaker ("greyer") dependence 
on wavelength than in a simple foreground screen model. 
Our model for dust absorption thus has essentially no sig- 
nificant free parameters. Results are sensitive mostly to the 
value of Tvo, which is calculated directly from our other 
model quantities. 

Our modelling of dust extinction is a major improve- 
ment over what has been done previously in semi-analytic 
models, both in terms of including a realistic 3D distribution 
for the stars and dust, and in terms of calculating the dust 
optical depth in a physically self-consistent way. The first 
sem i-anal ytic models to include dust were those of Lacey et 
al. (1993), using the dust and st ellar p opulation model of 
Guiderdoni & Rocca-Volmerange (1987), and Guiderdoni et 
al. (1998). They modelled the star and dust distributions 
as a uniform ID slab, but calculated the dust content self- 
consistently from a closed-box chemical evolution model. 



Kauffmann et al. ( 1999a ) and Somerville & Primack ( 199!: ) 



also use the ID slab model, but instead of predicting the 
slab optical depth, they use a power-law relation between 
dust optical depth and galaxy luminosity that is estimated 
from observations of z = galaxies. 

The main deficiencies of our current dust model are that 



it does not allow for clumping of the dust and stars or deal 
well with bursts, and that it only calculates absorption by 
dust, but n ot the spectrum of dust emission. However, Silva 
et al. (1998) have developed a more sophisticated dust model 
which includes both clumped and smooth components of the 
dust, deals accurately with bursts and is able to predict not 
only the extinction of starlight, but also the spectrum of the 
energy re-radiated by the d ust in the far infra-red and sub- 
mm. Granato et al. (2000) combine this dust model with 



our galaxy formation model to predict galaxy luminosity 
f unctio ns in the far infra-red and sub-mm and Lacey et al. 

investigate the high redshift behaviour and predict 



( ^ooo4 



number counts and integrated radiation backgrounds. 



5.4 Emission line Modelling 

We model the emission lines from photo-ionized gas in our 
galaxies by calculating the luminosity in Lyman continuum 
photons from the stellar population using the Bruzual & 
Chariot models, and combining this with HII region models 
to calculate line luminosities and equivalent widths. We have 
calculated results for important lines used as star formation 
indicators, such as Ha and OIL This is described in detail 



in a separate paper (Lacey et al. 2000b) 



6 METHODOLGY 
6.1 Model Parameters 

The complete hierarchical model of galaxy formation de- 
scribed in the previous section contains a significant number 
of parameters. However, relatively few should be considered 
as free parameters. These fall into three distinct categories: 
numerical parameters, parameters of the cosmological model 
and, finally, parameters related directly to our modelling of 
the physics of galaxy formation. 

The parameters that fall in the first set include the mass 
resolution, Mros, the number of timesteps in the merger tree, 
A'stcps and the starting redshift, Zstart. These do not rep- 
resent freedoms of the model, and we must simply choose 
values such that the quantities of interest have converged 
and are insensitive to further improvements. Also, there are 
options such as adopting singular isothermal spheres to de- 
scribe the dark matter and gas density profiles, or varying 
the distribution of halo spin parameters, which are not to 
be viewed as viable alternatives. Instead, we have included 
them simply in order to be able to vary our assumptions so 
as to gain insight into why the model behaves in a particular 
way. In these examples, any viable model should employ the 
options that are consistent with results of the high-resolution 
simulations that we are trying to emulate. 

The second set are parameters that specify the back- 
ground cosmological model. These include the density pa- 
rameter, ilo', the cosmological constant, Aq; the Hubble con- 
stant, h; the baryon density, Qh', and the shape and am- 
plitude of the linear theory mass power spectrum, P{k). In 
principle, each of these can be determined from observations 
that do not depend on galaxy properties. For example, most 
of these cosmological parameters are likely to be determined 
accurately from microwave background anisotropy measure- 
ments to be carried out by the MAP and Planck satellite 
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missions (e.g. Bond, Efstathiou & Tegmark 1997). Alter- 



natively, the power spectrum amplitude, ag, may be fixed 
by reference to the abundance of galaxy clusters, while the 
baryon density may be constrained by models of primordial 
nucleosynthesis and the observed abundance of the light ele- 
ments, like deuterium, at low and high redshift. Our general 
approach is to set these parameters according to such ex- 
ternal constraints. However, some properties of the galaxy 
formation models are particularly sensitive to the baryon 
density, Q,h, and to the normalization, erg. For this reason, we 
sometimes allow some variation of these parameters around 
the values otherwise indicated by the external observational 
constraints. 

The final set of parameters are those with which we di- 
rectly characterize our physical model of galaxy formation. 
Firstly, there is the IMF and its associated yield of metals, 
p, and the fraction, R, of stellar mass that is liberated in 
stellar winds and SNe. In principle, p and R are fixed by 
the choice of IMF, but, in practice, although R is quite well 
constrained, p is quite uncertain. Secondly, we have the pa- 
rameters, e*, a-t, Vhot and ahot in the star formation and 
feedback laws. Then, there is the parameter, e, the frac- 
tion of the metals produced in SNe which escape directly to 
the hot diffuse halo gas. Also, we require the vertical scale- 
height of dust relative to that of the disk stars (although 
our results are very insensitive to this parameter). Finally, 
there are the parameters /df and /eiup, which modulate the 
frequency of galaxy-galaxy mergers and determine when a 
merger results in the formation of a spheroid. Although the 
number of model parameters is not small, we shall see that 
the resulting freedoms of the model are still quite limited 
and that only a small subset of observed properties of low 
redshift galaxies are needed to constrain a model fully. 

6.2 Model Output 

The output of our code is a list of the galaxies that form in 
each simulated halo at one or more redshifts. For each galaxy 
the output lists: a flag which indicates whether the galaxy 
is the central galaxy of the halo in which it is contained or 
a satellite galaxy; the mass of cold gas in the disk; the mass 
of stars in the disk and bulge; the luminosities in any chosen 
band of the stars in the disk and bulge; selected emission line 
luminosities and equivalent widths; the half-mass radius of 
the disk and bulge individually and combined; the circular 
velocities at the half-mass radii of the disk and bulge; the 
metallicity of the cold gas and also, if the galaxy is a central 
galaxy, the metallicity of its hot gas halo; the metallicities 
and age of the bulge and disk stars weighted by mass or 
by luminosity in any selected band; the instantaneous star 
formation rate in the disk; the mass and circular velocity at 
the virial radius of both the halo in which the galaxy was 
last a central galaxy and the halo in which it is contained 
at the chosen output redshift. The effect of dust within each 
galaxy on the luminosities and line strengths is computed 
in an additional step by assuming an inclination angle for 
each galaxy. Also, since we know the disk and bulge sizes, we 
can compute surface brightness distributions and isophotal 
magnitudes for each individual galaxy, assuming exponen- 
tial profiles for disks and r^^* profiles for spheroids. Since we 
also know the number density of each of the halos we have 
simulated, it is straightforward to estimate galaxy luminos- 



ity functions and galaxy number counts (both using either 
total or isophotal magnitudes) or to sample our output to 
build up either volume-limited or magnitude-limited galaxy 
catalogues from which to draw galaxy samples for compari- 
son with observational datasets. In this work we have used 
the halo abundance given by the Press-Schechter formal- 
ism, but it is now possible to adopt improved analytic esti- 
mates (Sheth, Mo & Tormen ^OOOj ) which accurately match 
the abundance of ha los found in large N-body simulations 
(Jenkins et al. 200C). We have checked that switching to 
these more accurate formulae changes the model results far 
less than varying some of the galaxy formation parameters. 

Once we have calculated a model, we have the ability 
to select from the output any particular galaxy and recom- 
pute its formation history, this time choosing to output its 
properties more frequently and to record its star formation 
and merger history. In this way we can generate the com- 
plete formation history of selected galaxies and, for each, 
construct their own individual merger trees. Examples of 
galaxy merger trees constructed in this w ay have been pre- 
sented in Figure 9 of Baugh et al. (1998). 



6.3 Strategy 

Our adopted methodology is to select a cosmological model 
based on constraints from large-scale structure and then 
vary the galaxy formation parameters in order to match 
as best as possible a selection of low-redshift observational 
data. Since galaxy formation is undoubtedly a complex pro- 
cess, the simple model we have constructed cannot aspire 
to be a complete and full description. Thus, it is inevitable 
that in some cases our models will only produce a moderate 
level of agreement with certain observational data. 

In the following section, we illustrate this process of 
defining the parameters of the galaxy formation model for 
one particular cosmology. We use this example to illustrate 
the way in which we apply the observational constraints and 
to show how the model predictions depend on each of the 
model parameters. 



T OBSERVATIONAL CONSTRAINTS ON THE 
MODEL AND EFFECTS OF VARYING THE 
PARAMETERS 

In the following sub-sections, we compare models con- 
structed within a particular cosmology with a variety of 
statistics estimated from the observed properties of the local 
galaxy population. For each statistic, we illustrate how the 
predictions of the galaxy formation model depend on the pa- 
rameters, and present one model, our reference model, which 
is the best compromise when measured against a full range 
of observational data. The observational constraints used to 
fix each of the main parameters of our reference model are 
as follows: 

• Qfhot : faint end of luminosity function and Tully- Fisher 
relation 

• Vt^ot- faint end of luminosity function and sizes of low 
luminosity spirals 

• £*: gas fraction for spirals 

• a*: variation of gas fraction with luminosity 
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• /ciiip: morphological mix for galaxies 

• IMF: observations of solar neighbourhood 

• T: in luminosity function 

• p: metallicity of L* ellipticals 

The chosen parameters values of our reference model are 



listed itl iable i. 



We emphasize that not all of the observational data pre- 



sented iln this section are used to tlx model parameters - some 
of the data fttovide tiistS 6f th6 m6d«!l, iind aid ah6wh ih thiri 
section to illustrate the effects of varying the parameters. 

The cosmology that we have chosen in order to illustrate 
how observational data may be used to constrain the galaxy 
formation parameters is a flat, low-density cold dark matter 
model with a cosmological constant. The parameters that 
we adopt for this ACDM model are Oq — 0.3 and Aq = 0.7. 
Such a model is currently favoured by quite a range of ob- 
servational evidence. For reasonable values of the Hubble 
constant {h ~ 0.7), the shape of the mass power spectrum 
is in good agre ement with estimate s from large-scale galaxy 
clustering (e.g. Maddox et al. 199C). The value of Qq is con- 



sistent w ith the high barvon fr a ction in clusters (White et 
al. 199| [White fc Fabian 1995|; Mohr & Evrard 1997|) and 



with the mild ev olution in the abundance of X-ray clusters 
(Eke et al. 1998). The joint values of flo and Ap are in ac- 
cord with estimates from h i gh redshift SNe (Perlmutter et 
al. 199|[ [Riess et al. 1998| ; [Garnavich et al. 199^ ), while 
f2o + Ao = 1 is in ag reement with the detection of the first 
CM B floppier peak ([de Bernardis et al. 200C ; Lange et al 
200C |; |H|annay et al. 2000^ [Balbi et al. 200C| ). For this model 



the normalization of the mass power spectrum, erg, derived 
from the number density of X-ray emitting galaxy clusters 
is consisten t with that from the amplitude of CMB fluctu- 
ations (e.g. Cole et al. 1997), and both are consistent with 
the po wer sp ectrum of the mass a t z — 2.5, inferred by Croft 
et al. (19991) and Weinberg et al. (1999) from analysis of the 
Lyman- a forest. 

The galaxy formation model as a whole is a complex 
system, with the result that the dependence of a particular 
statistic on a given parameter can be complicated. Thus, 
when one parameter is varied, the behaviour of the model 
certainly depends on the other constraints that have been 
applied. This fact, combined with the number of parame- 
ters, makes it unfeasible to present the full range of possible 
model behaviour. One should therefore be careful not to 
over-interpret the trends discussed below: since the effects 
of varying one parameter can depend on the values of the 
others, it is dangerous to assume that these trends can be 
used to assess the result of varying more than one parameter 
at a time. 

As well as varying the parameters that regulate the 
physics of galaxy formation, we allow some variation of 
the parameters, as, h and fib, that define the cosmological 
model. However, these are only allowed to vary within the 
ranges permitted for consistency with estimates of the abun- 
dance of rich galaxy clusters, Hubble's constant, and pri- 
mordial nucleosynthesis. Whenever we vary any of these pa- 
rameters, we consistently adjust the power spectrum shape 
parameter, F, accord ing to the fitting formula for CDM pro- 
posed by Sugiyama (1995): 



The estimated values and 1-a errors that we a dopt for these 
Quantities are as = 0.93 d z0.07 (|Eke et al. 1996|) , h = 0.7±0.1 
( Preedman et al. 1998[ [Madore et al. IQQSj ), and J^b^^ = 
0.0125±0.0025 ( [Walker et al. 1991| ). We note that somewhat 
higher values of Qhh are favoured by recent est imates of 
the D/H ratio f rom QSO absorption line systems ( ^chramm 



Turner 1998| ; [Buries fc Tytler 1998[ ), and b y models of 
the mean optical depth of the Lyman-a forest ( [Ranch et al 



1997[ ; [Weinberg et al. 1997| ). We discuss the consequences of 



increasing our adopted value of S7b in Section ^. 

In a similar fashion, many of the parameters that de- 
scribe the physics of galaxy formation can only plausibly 
be varied by modest amounts. Consider, for example, the 
threshold, /oiiip, above which galaxy mergers are termed vi- 
olent and assumed to result in the formation of an elliptical 
galaxy. By definition, fdUp < 1, and based on simulations 
of galaxy me rgers, a reas onable lower limit is /dup <; 0.3 
(Walker et al. 1996; Barnes 1998). Other parameters that are 
similarly constrained to lie within relatively narrow ranges 
are fdf and, to a lesser extent, p. The merger timescale co- 
efficient, /df, should be close to unity if an infalling galaxy 
retains its dark matter halo throughout most of the time 
when dynamical friction is removing angular momentum and 
energy from its orbit. It could be larger than unity if the 
dark matter halo is efficiently stripped off at early times, 
but cannot plausibl y be significantly smaller than unity. As 
described in Section 



5.1 



the choice of IMF quite accurately 
determines the recycled fraction, R, and also sets some con- 
straint on the yield, p. 

7.1 Galaxy Luminosity Functions 

The single most important constraint on our galaxy forma- 
tion models is the local galaxy luminosity function. This is 
one of the most fundamental properties of the galaxy popu- 
lation and it is also one of the best measured, at least over 
a restricted range of luminosities. Matching the galaxy lu- 
minosity function is a prerequisite for any realistic model of 
galaxy formation, because the detectability of galaxies de- 
pends directly on their luminosity. Thus, a model that fails 
to match the bright end of the luminosity function can lead 
to misleading conclusions when tested, for example, against 
samples selected by apparent magnitude. 

Fig. ^ shows estimates of the local galaxy luminosity 
function in the blue optical bj-band and in the near infrared 
K-band. The bj-band data are taken from the APM-Stromlo 



F = no/iexp(^-^(V2ft-|-!^o)) • 



(7.1) 



galaxy r edshift survey (Loveday et al. 1992 ), th e ESQ slice 
project (Zucca et al. 1997), the DUKST survey (Ratcliffe et 
al. 1998 ) and fr om preliminary resul ts from the 2dF galaxy 
redshift survey ( Maddmc et al. 1998 ). The K-ban d dat a are 
from Mobasher et al. (199S), Glazebroook et al. (1995) and 
Gardner et al. (1997). In both bands, there is reasonably 
good agreement among the various estimates at the bright 
end. At the faint end, however, there is considerable disper- 
sion among the results from different surveys. In the bj-band 
these differences are large compared to the statistical errors. 
We must, therefore, conclude that either the galaxy luminos- 
ity function differs in the different volumes surveyed, or that 
systematic differences in the selection characteristics of the 
surveys give rise to these differences. 

The solid and dashed model curves shown in Fig. ^ cor- 
respond to the reference model, both with and without the 
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Table 1. The values of the model parameters for the reference model. 

Ao f^b r o-g a* Vhot «hot e /dlip /df IMF p _R T 

0.3 0.7 0.02 0.7 0.19 0.93 0.005 -1.5 200.0 2.0 0.0 0.3 1.0 Kennicutt 0.02 0.31 1.38 




Figure 5. Comparison of the bj and K-band galaxy luminosity functions in the ACDM reference model with a compilation of obser- 
vational data. The solid line is the model including the effects of dust, with T chosen so as to obtain agreement with the observed 
luminosity functions at M^j — 51og/i = —19.8. The dashed lines show the corresponding luminosity functions before the effects of dust 
extinction are included. The dotted line on the left hand panel shows how the luminosit y fu nction is modified if the isophotal magnitude 
within an isophote of 25 mag/arcsec^ is used instead of the true total magnitude. See §[7.3| for details. 



inclusion of dust. In the model with dust, the vertical scale- 
height of the dust distribution was taken to be the same as 
that of the disk stars, but varying this assumption makes 
very little difference to the luminosity functions. The model 
without dust produces a reasonable K-band luminosity func- 
tion, but it gives rise to galaxies which are systematically too 
bright in bj . By contrast, the model with dust is a good fit to 
the bright end of both the bj and K-band luminosity func- 
tions. We shall see below that the differential effect of dust, 
which generates greater extinction at shorter wavelengths, is 
important in providing a good match to the observed B~K 
colour distributions. It also results in a model that comes 
significantly closer than our previous models to matching 
simultaneously the zero-point of the observed I-band TuUy- 
Fisher relation and the bright end of the bj-band luminosity 
function, thus largely overcoming an im portant shortcoming 
of our earlier models (Cole et al. 1994). The importance of 



the role of dust in achieving this si multan eous match has 
also been no ted by Kauffmann et al. ( 1999a ) and Somerville 
& Primack (1999). The prescription for the dust distribu- 



tion assumed in Fig. |^ will be retained in all the following 
comparisons. 

Many of the model parameters have a direct effect on 
the galaxy luminosity function. Typically, varying any one 
of them on its own causes only a minor change in the shape 
of the luminosity function, but can cause a significant over- 
all shift (either left or right) in the luminosity scale. Our 



normalization strategy separates out these two effects by 
adjusting the parameter T so as to keep the amplitude of 
all the luminosity functions fixed at Mb, = —19.8 + 5 log ft. 
Recall that is the initial stellar mass fraction in lumi- 



nous stars (equation 5.2), and that the remaining fraction 
is assumed to be made up of non-luminous brown dwarfs. 
Physically, one ought to vary the net recycled fraction R 
when adjusting T, so as to keep the value 7?i = TR constant 
for the luminous stars alone (§ ^.2| ). The reference model has 
T and R chosen consistently to give _R i — TR — 0.42 for 
the Kennicutt IMF, as found in Section 5.2, To achieve this 
requires some iteration as the value of T is determined after 
the model has been run by matching a point in the B-band 
luminosity function. For this reason, when showing the ef- 
fects of varying parameters we choose to keep R constant 
rather than Ri. 

The effect of parameter changes on the position of the 
luminosity function is summarized in Table 2. Here, we give 
the values of T required to shift each lu minos ity function 
into coincidence with the Zucca et al. (1997) luminosity 
function at Mbj = -19.8-1-5 log h. We will see in Section 7.7 
that in models which have realistic stellar mass-to-light ra- 
tios, T ~ 1.3 — 2 for the Kennicutt IMF. In a few cases, how- 
ever, this normalization procedure requires T < 1, which 
is unphysical. (It corresponds approximately to removing 
from the IMF all of the brown dwarfs and some of the low 
mass visible stars, by increasing the lower mass limit above 
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Figure 6. The effect on the bj-band luminosity function of vary- 
ing the star formation and feedback parameters, «hot a-^d ^hoti 
and the assumed halo DM and hot gas density profiles. In all 
cases, the value of T has been fixed by requiring the model lumi- 
nosity functions to agree with the observations at j — 5 log h = 
— 19.8. a) Shows how increasing Ohot suppresses the formation 
of low-luminosity galaxies and so controls the faint-end slope 
of the luminosity function. The dotted curve is for Qhot = Ij 
solid for Ohot = 2, and dashed for Qhot = 5.5. b) Demonstrates 
how increasing Vhot lowers the faint end of the luminosity func- 
tion. Results are shown for Vhot = lOOkms"^ (dotted curve), 
^hot = 200kms~^ (solid curve), and Vhot = 300kms~^ (dashed 
curve), c) Shows how the bright end of luminosity function de- 
pends on the model adopted for the density profile of the hot gas. 
The solid curve is our reference model which has an NFW profile 
for the DM and the "/3-model" for the gas, with a core radius that 
dep ends o n the fraction of gas that has previously cooled (see Sec- 
tion ^!l^). The dotted and long-dashed curves also assume NFW 
profiles for the DM, but assume a fixed core radius "/3-model" 
(dotted) or an NFW profile for the gas (long-dashed). The short- 
dashed line is for a singular isothermal sphere (p oc r~^) model 
for both gas and DM. 



O.IAIq.) Nevertheless, we present these models because they 
help to clarify some of the trends that occur when a single 
parameter is varied. For the derived normalization of T, the 
table gives other properties of the galaxy population. Sev- 
eral of the entries in this table are discussed further in the 
relevant sections below. 

Two parameters that strongly affect the shape of the 
galaxy luminosity function are Vhot and cthot, the qua ntitie s 
that define our model of stellar feedback (equation 1.15 ). 
As can be seen in Fig. the faint-end slope of the lumi- 
nosity function is very sensitive to ahot, with large values 
producing a shallower slope. Similarly, as Fig. ^ shows, in- 
creasing Vhot also reduces the number of faint galaxies. Both 
these dependencies are easily understood: stronger stellar 
feedback makes it increasingly more difficult for luminous 
stars to form in low-mass halos. To match the ver y shal low 
faint-end slope s een i n the data of Loveday et al. ( 1992 ) or 
Ratcliffe et al. (1996), requires a hig h value of ahot, such 
as that adopted by Cole et al. (199'!l), who compared their 
models against the first of these surveys. Such a value, how- 
ever, would lead to a disagreement with the data of Zucca 
et al. (1997). The differing observational estimates of the 



luminosity function indicate that the faint-end slope is not 
as robustly determined as one might wish, perhaps because 
it depends on the details of the survey selection criteria. We 
therefore do not use it as a model constraint. Instead, we will 
see in Section that extreme values of a^ot are disfavoured 
on other grounds and this leads us to favour models whose 
luminosity functions have quite steep faint end slopes. 

The bright end of the luminosity function is sensitive 
to the density profile assumed for the halo gas, because this 
controls how much of the gas can cool. This effect is not 
important in low-mass halos, in which the gas temperature 
Tvir is low enough that most of the gas can cool anyway, but 
it becomes important in large groups and clusters, in which 
only the dense central regions have time to cool. Fig. ^c com- 
pares the effects of using different gas profiles. Our reference 
model (shown by the solid line) assumes an NFW dark halo 
and a /3-model for the gas (equation 4.2), with a core radius 
that starts at rcoro = ^new/S and grows depending on how 
much gas has already cooled in progenitor halos. The model 
luminosity function and other properties are not sensitive to 
the precise value of this initial core radius. For example, if 
instead, we set Tcorc = ^new/O as the initial value, then the 
change in the luminosity function is almost too small to be 
visible and the other properties listed in Table 2 also only 
vary slightly. In principle, constraints can be placed on the 
initial gas density profile from the observed X-ray emission 
profiles of groups and clusters, but in practice this requires 
complex modelling to take account of the emission associ- 
ated with the central cooling flow. 

Our model fits the observed bright end of the luminos- 
ity function well. It is compared in the figure to a model in 
which the gas core radius is kept fixed at rcorc = '"nfw/3 
(dotted curve), another in which both gas and dark matter 
have the same NFW profile (long-dashed curve), and finally 
to a model in which both gas and dark matter have singu- 
lar isothermal sphere profiles (short-dashed curve), as has 
been assumed in most previous work. The latter three mod- 
els produce many more high luminosity L ^ galaxies 
than is observed. This difference in the assumed halo gas 
profiles explains most of the differences in the shape of the 
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Table 2. The variation of the average colour, the mass-to-hght ratio of the stellar populations and the zero-point of the TuUy-Fisher 
relation with model parameters. In all cases, T is adjusted so as to match the observed bj-band luminosity function at -A/bj = —19.8 + 
5 log h. The first column lists the parameter that has been varied relative to the reference ACDM model. The second column gives the 
required value of T. The third lists the median B-K colour for galaxies in the range —24.5 < Mk — 51og/i < —23.5. The following 
four columns list the median B- and I- band stellar mass- to-light ratios, in units of ^Ma/L©, for disk and elliptical galaxies with 
—20 < Mb — 51og/i < —19.0. These are compared with observed values in Section 7.7 The last two columns show the offset in 



magnitudes from the observed I-band TuUy-Fisher relation at Vdisk = 160 km s , and the median ratio of the circular velocity of the 
disk to that at the virial radius of the halo in which it formed for galaxies with —20 < Mi — 51og/i < —18. 



Modified Parameter 


T 


B-K 


Disk: M/L-B 


Disk: M/Li 


Elliptical: M/L-g, 


Elliptical: M/Lj 


AM/ 


v A ; r^!.. / vv. ^ 1 ^ 
* clisK / * nalo 


Reference Model 


1.38 


3.86 


2.0 


1.8 


5.4 


2.9 


0.98 


1.35 


No Dust 


2.30 


3.60 


2.1 


2.4 


8.7 


4.9 


1.48 


1.41 


= 1 


1.39 


3.81 


2.1 


1.8 


4.8 


2.8 


1.13 


1.50 


Qvifit =5.5 

not ■ 


1.1 


4.01 


1.4 


1.2 


4.3 


2.3 


1.02 


1.02 


Vhot = lOOkms"^ 


1.32 


4.11 


2.5 


2.0 


6.3 


3.1 


1.73 


1.71 


Vhot = 300kms-i 


1.15 


3.56 


1.2 


1.2 


2.5 


1.8 


0.75 


1.20 


a* = 0.0 


1.28 


3.96 


2.0 


1.7 


5.4 


2.8 


0.83 


1.27 


a* = -2.5 


1.31 


3.81 


1.7 


1.5 


5.0 


2.8 


1.09 


1.41 


rJb = 0.01 


0.45 


3.62 


0.5 


0.5 


1.0 


0.7 


0.23 


1.13 


= 0.04 


3.07 


4.11 


7.3 


4.7 


13.9 


6.9 


2.12 


1.96 


h = 0.6 


0.95 


3.87 


1.7 


1.5 


4.8 


2.6 


0.92 


1.34 


h = 0.8 


1.77 


3.86 


2.1 


1.9 


5.4 


3.0 


1.19 


1.38 


e* = 0.01 


1.70 


3.85 


2.1 


1.9 


7.4 


3.9 


1.12 


1.30 


e* = 0.0033 


1.13 


3.90 


1.8 


1.5 


4.5 


2.5 


1.03 


1.40 


p = 0.0075 


1.66 


3.28 


1.6 


1.7 


3.9 


2.7 


1.01 


1.36 


p = 0.03 


1.17 


4.16 


2.1 


1.7 


5.7 


2.8 


0.92 


1.34 


R = 0.19 


1.31 


3.80 


2.1 


1.9 


5.9 


3.2 


1.07 


1.38 


R = 0.49 


1.41 


3.97 


1.7 


1.4 


3.7 


2.0 


0.85 


1.31 


as = 0.86 


1.43 


3.82 


2.0 


1.7 


4.2 


2.5 


1.03 


1.34 


as = 1.0 


1.48 


3.90 


2.2 


1.9 


6.1 


3.3 


1.20 


1.38 


IMF: Salpeter, H = 0.28 


0.79 


3.85 


1.9 


1.7 


5.5 


3.0 


1.00 


1.35 


/form = 1.5 


1.27 


3.82 


1.9 


1.6 


5.1 


2.8 


0.93 


1.42 


/df = 0.5 


1.34 


3.89 


2.0 


1.7 


4.4 


2.5 


1.03 


1.36 


/df = 2.0 


1.16 


3.80 


1.5 


1.4 


4.8 


2.6 


0.84 


1.34 


/ellip =0.5 


1.38 


3.90 


2.1 


1.8 


5.8 


3.1 


0.99 


1.35 


^core = rNFw/6 


1.40 


3.85 


2.0 


1.8 


5.4 


3.0 


1.05 


1.35 


Fixed gas core radius 


1.23 


3.97 


2.3 


1.8 


2.9 


1.9 


0.87 


1.34 


NFW gas traces DM 


1.23 


4.04 


2.5 


1.8 


3.1 


2.0 


0.92 


1.34 


SIS gas traces DM 


1.08 


4.24 


2.4 


1.7 


3.4 


2.0 


0.95 


1.41 


Unstable disks 


1.50 


3.91 


2.1 


1.9 


4.7 


2.8 


1.11 


1.28 


Accretion by disk 


1.40 


3.89 


2.3 


1.9 


5.5 


3.0 


1.12 


1.38 



bright end of the luminosity function betw een our refer ence 
model and t he mo dels of Kauffmann et al. (1993,1999a) and 



As was shown in Cole et al. (1994), the shape of the 



Somerville (1997), although the procedure in these papers 



of using the Tully-Fisher relation rather than the luminosity 
function as the primary observational constraint also has an 
effect. These authors all invoke an artificial cut-off on the 
circular velocity of halos in which gas is allowed to cool to 
form visible stars, in order to ameliorate their problems in 
fitting the luminosity function. We believe that our standard 
model for cooling is the most physical ly rea sonable in this 
regard, for the reasons given in Section 4.1.1, More recently, 
Somerville & Primack ( |l99c| ) have presented models with no 
artificial cooling cut-off which nevertheless produce a good 
match to the bright-end of the B-band luminosity function. 
In this case, this improvement is achieved partly as a result 
of the empirical dust model they have adopted, which has 
an extinction that increases with increasing galaxy luminos- 
ity. However, as dust has less effect in the K-band, they find 
that for some of their models, the shape of the bright-end 
of the K-band luminosity function remains a poor match to 
observations. 



luminosity function is also influenced by the efficiency of 
galaxy mergers, which is controlled by /df. However, we now 
impose the const raint /df <: 1, a limit suggested by numer- 
ical simulations (Navarro et al. 1995a), which also leads to 



an acceptable morphological mix in the model. With this 
bound, the residual variation in the shape of the luminosity 
function with /df is small compared to its dependence on 
the feedback parameters Viiot and tthot and on the halo gas 
profile. Our treatment of feedback is, in fact, the main fac- 
tor responsible for the overall shape of the model luminosity 
function at L L*, and our assumptions about cooling are 
the main determinant of the shape at L ^ L*. 



7.2 The Tully-Fisher Relation 

In Fig. ^, we compare our predicted I-band Tully-Fisher 
(TF) relation with the observed relation defined by the com- 
plete diameter-limi ted s ubset of spiral galaxies selected by 
de Jon g fc Lacey (2000) from the catalogue of Mathewson 
et al. (h99a). The observed circular velocity plotted here is 
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Figure 7. The dependence of the model I-band TuUy-Fisher re- 
lation on a) the sample selection criteria, b) the feedback pa- 
rameter, cihoti S'lid c) the baryon density, Cj,. In each case, the 
model curves trace the median magnitude as a function of cir- 
cular velocity, and the errorbars show the 10 and 90 percentiles 
of the distribution. The magnitudes are face-on values, includ- 
ing the effects of dust. The points show the observed distribution 
for a subsample of Sb-Sd galaxies selected by d e Jon g & Lacey 
(2000) from the Mathewson, Ford & Buchhorn( 1992 ) catalogue 



and, again, all magnitudes have been corrected to face-on values. 
All the curves in the top panel are for the reference model, but for 
different galaxy selection criteria. The dotted line is for all spiral 
galaxies which are the central galaxies in their halos, the dashed 
Une for all spiral galaxies (both central and satellite), and the 
solid line for all star-forming spiral galaxies with gas fractions of 
10% or greater (this final selection is retained in b & c). The thick 
solid line shows, for central galaxies, the result of using the cir- 
cular velocity at the virial radius of the halo in which the galaxy 
formed, rather than the disk circular velocity. In b), the solid line 
refers to the reference model (which has cthot = 2), while the 
dashed line is for Ohot = 1 ^^nd the dotted line for a^ot = 5-5- In 
c), the solid line refers, again, to the reference model (which has 
f2b = 0.02), while the dashed line is for Q^, = 0.01 and the dotted 
line for Ob = 0.04. 
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the maximum, Vmax, of the measured rotation curve. The 
observed I-band magnitudes have been corrected to face- 
on values. For the model, we use the dust-extincted I-band 
magnitude for galaxies seen face-on and the circular veloc- 
ity, Vdiak, evaluated at the half-mass radius of the disk. Vdiak 
includes the self-gravity of the disk, and is evaluated in the 
disk mid-plane, as discussed in Appendix |^. The peak of 
the rotation curve may occur at a radius other than the 
half- mass radius, in which case the quantity, Vdisk, that we 
plot may be systematically low compared to the measured 
Vmax- However, we expect the difference to be small since 
our model galaxies typically have reasonably flat rotation 
curves, as indicated by the values of Vdisk/ Vhaio listed in 
Table 2, which are close to unity. 

The upper-panel in Fig. ^ shows how the model TF 
relation depends on the sample selection criteria. In all 
cases, we have selected disk-dominated galaxies with dust- 
extincted I-band bulge-to-total light ratios in the range 0.02 
to 0.24, to match approximately the range of ga laxy types, 
Sb-Sd, selected in the Mathewson et al. (1992) catalogue. 



This makes use of the approximate conversion between Hub- 
ble T-type, (T = 3-7 for Sb-Sd) and bulge-to-disk ratio, 
described in Baugh et al. ( 1996b ) and based on the data 
of Simien & de Vaucoleurs (1986). The model TF relation 
for this cornplete galaxy sample is shown by the dashed 
line in Fig. Ma. The slope of the predicted TF relation is 
close to that of the data, and this remains true for the 
sub-samples discussed below. It has an offset of 1.2 mag- 
nitudes relative to the zero-point of the observed relation 
and a spread between the 10 and 90 percentiles of the dis- 
tribution at Vc — 160 kms^^ of 1.7 magnitudes. This spread 
is significantly larger than that for the observational sample, 
which is 1.1 magnitudes. 



In Cole et al. ( |l994|) the TF relation was plotted for 
galaxies at the centres of halos only. Here the result of se- 
lecting only central galaxies is shown by the dotted line. 
The exclusion of satellite galaxies removes some galaxies 
which have exhausted their reservoirs of cold gas and so 
have faded as their stellar populations have aged. This has 
the effect of producing a somewhat tighter TF correlation, 
with a spread of only 1.2 magnitude, and of reducing the off- 
set in the zero point to 1.0 magnitudes (at Vc = 160 kms~^). 
A more realistic selection is to consider only disk galaxies 
with a significant cold gas fraction, which we take to be 
Mgas/(Mgas -I- Afstars) > 10%. This Is reasonable, as without 
ongoing star formation, disk galaxies will not have promi- 
nent, recognizable spiral arm features. In addition, inter- 
stellar gas is required for the measurement of the rotation 
velocity in TF datasets, either to produce the emission lines 
from which optical rotation curves are measured, or to pro- 
duce the HI emission used in HI rotation measurements. The 
TF relation for this subsample of star-forming spiral galax- 
ies is shown by the solid curve in Fig. and is repeated as 
the solid curve in the lower two panels. It has a spread of 
0.98 magnitudes, which is slightly smaller than the observed 
spread of 1.1 magnitudes. The offset in the zero point of the 
relation is 0.98 magnitudes, which is equivalent to a factor 
of approximately 1.3 in circular velocity. Since the effective 
mass-to-light ratio in our models is normalized (through T) 
by reference to the bright end of the bj-band luminosity 
function, we find that both the zero-point and the scatter 
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in the model TuUy-Fisher relation are insensitive to most 
changes in the galaxy formation model parameters. 

The parameters that do have an effect are Ohot and 
S7b- This is illustrated by the curves in Fig. ^ and 7c re- 
spectively. Increasing the feedback parameter, Qhot, makes 
it increasingly difficult to form stars in low-circular velocity 
galaxies. Consequently, the luminosity of low-Vdisk galaxies 
is reduced, and the model TuUy-Fisher relation bends away 
from the observed correlation at faint magnitudes. A value of 
Qhot ~ 2 is required to produce a correlation that runs paral- 
lel to the observed relation over the full range of magnitudes 
probed by the data. The effect of increasing fib (Fig. ^) is 
to cause the model TuUy-Fisher relation to bend away from 
the observations at bright magnitudes. The reason for this is 
that, in order to maintain a match to the bright end of the 
b,j-band luminosity function, larger fib requires a larger T 
and thus larger mass-to-light ratios for all the galaxies. The 
self-gravity of bright spiral disks then plays a larger role in 
determining the galaxy's rotation curve. This is quantified 
by the ratio, Vdisk/Vhaio, listed in Table 2, which increases 
substantially as fib is increased. It is this effect that leads 
us to favour a relatively low value of Q,h as compared to the 
currently most favoured numbers derived from primordial 
nucleosynthesis considerations. Note that a very low value, 
f2b = 0.01, results in a good match to the zero-point of 
the TuUy-Fisher relation. However, this apparent success is 
at the cost of an unphysical value, T = 0.49, required to 
make galaxies bright enough to match the bj-band luminos- 
ity function. 

The failure of our model to produce a TuUy-Fisher re- 
lation with a zero-point that matches the observations well 
is reminiscent of a similar short coming in the earlier Q. = 1 
CDM model of Cole et al. ( 1994 ) and also, b ut to a lesser ex- 



tent, the low-fio CDM models of Heyl et al. ( 1995 ) . However, 
this discrepancy hides a significant improvement in our new 
models. In our previous work we did not attempt to model 
the internal mass distribution within a galaxy, and simply 
took the circular velocity of the galaxy to be that at the 
virial radius in the halo in which it formed. If we followed this 
same procedure now, and plotted V^^io rather than Vd isk as a 
function of I-band magnitude, we would find a near-perfect 
match to the observed TuUy-Fisher relation, as indicated by 
the heavy solid line in Fig. Ma. The reason for this difference 
in the Vhaio-A^/ relation between our old and new models is 
largely the inclusion of dust in the new models, which helps 
in two different ways to simultaneously match the bj-band 
luminosity function and the I-band TuUy-Fisher relation. 
Firstly, dust makes galaxies dimmer in b,j, allowing a better 
match to the observed luminosity function with a smaller 
value of T, but it also makes them redder in bj-J. The net 
effect is that the I-band luminosities used in the TF relation 
are increased. Secondly, dust affects the calculation of the 
luminosity function and the TuUy-Fisher relation in differ- 
ent ways, because observational estimates of the luminosity 
function use magnitudes uncorrected for dust, whereas ob- 
servational estimates of the TuUy-Fisher relation partially 
correct for the effects of dust through the correction to face- 
on magnitudes. Some of the se effects are also discussed by 
Somerville & Primack (1999). Increasing the amount of dust 
beyond that present in our reference model by, for instance, 
increasing the assumed yield, can further improve the TuUy- 



Fisher zero-point. However, this is achieved at the expense 
of making the galaxy population too red. 

7.3 Disk Sizes 

In our model, the sizes of galaxy disks are fundamentally 
determined by the angular momentum gained by the proto- 
disk material through the action of tidal torques, which are 
most effective when halos are turning around and collaps- 
ing, prior to becoming virialized. The distribution of halo 
spin parameters is well understood, and is reasonably ac- 
curately modelled by the distribution (equation 3.7) that 
we have adopted. The main sources of uncertainty are the 
distribution of angular momentum within a halo, which de- 
termines the angular momentum of that fraction of the gas 
that cools to form a disk, and whether the gas conserves its 
angular momentum during the coll apse. Th e assu mptions 



that we have discussed in Sections 3.2.3 & 4.1.3 



sonable, but are not direc tly su pported by simulations (see 



the discussion in Section 4.1.3) and warrant further inves- 
tigation. Apart from this, the largest remaining influence 
on the distribution of galaxy disk sizes is the strength of 
stellar feedback. If feedback is weak, stars form efficiently 
in small, dense halos at high redshift, while if feedback is 
strong, star formation is suppressed until larger halos form 
at lower redshift. Thus, increasing the value of Vhot results 
in galaxies having larger disk scalelengths at a given lumi- 
nosity. This dependence is shown explicitly in Fig. |^ and 
is weak for L ^ L*, but is significant at lower luminosities. 
A value of Vhot = 200 km s^^ produces a model for which 
the position of the peak in the disk scalelength distribution 
of spiral galaxies at different luminosities is close to what 
is found observationally by de Jong & Lacey ( 2000 ). More- 
over, the predicted width of the distribution is quite similar 
to that observed, although somewhat broader. Our model 
does not predict a large population of bright galaxies with 
either extremely large or small scalelengths. 

The long-dashed line in Fig. |^ shows how the size dis- 
tribution of disks in the reference model could be modified 
by the effects of disk instability. In this var iant, we have 
checked the disk stability criterion, equation (4.1J), at each 
timestep and have taken the material from unstable disks 
with em < 1 and added it to th e sphe roid or bulge com- 
ponent. Like Mo, Mao & White ( |l998a| ) , we find that this 
depletes the small disk scalelength side of the distribution 
and produces a slightly better match to the observed distri- 
bution for L ~ L* disks. 

It is interesting to examine the effects of surface bright- 
ness selection on estimates of the galaxy luminosity func- 
tion. For galaxies in the reference model, we computed the 
difference between the total magnitude and the magnitude 
within an isophote of 25 mag/arcsec^ in bj, for a galaxy 
seen face-on, assuming that the dust attenuation factors are 
constant with radius for the bulge and disk components. We 
then assumed that this aperture correction is independent 
of the actual inclination angle of the galaxy, and estimated 
the resulting galaxy luminosity function. Because of these 
approximations, and because in real galaxy surveys some 
attempt is made to extrapolate to total magnitudes, the re- 
sulting model luminosity function cannot be quantitatively 
compared to observations. Nevertheless, the difference be- 
tween this estimate, shown by the dotted curve in Fig. ^ 
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Figure 8. A comparison of predicted spiral galaxy disk sizes with observations. The two panels show the distribution of disk exponential 
scalelengths /io (number density as a function of scalelength and absolute magnitude) in luminosity bins on either side of L*. The points 
with error bars and t he tr iangles are, respectively, observational data and 95% confidence upper limits for Sb-Sd galaxies from the work 
of de Jong & Lacey ( 200C ) , allowing for the dependence of the observational selection function on galaxy size and luminosity. The lines 
are the model results for varying Vhoti fof galaxies with {B/T)j < 0.24. The dotted, solid and short-dashed curves are for Vhot = 100, 200 
and 300 kms~^ respectively. The long-dashed curve is for a variant of the reference model in which unstable disks have been converted 
to bulges, as described in the text. 



and that based on the true total magnitudes gives an indi- 
cation of the sort of systematic errors that could be present 
in real surveys. The main effect of using isophotal magni- 
tudes is to produce a small faintward shift at the bright end 
of the luminosity function, and a larger change in the faint- 
end slope. The dependence on the surface brightness limit 
suggests that the faint-end slope derived from surveys se- 
lected from photog raphic plates could be artificially shallow 
(see also McGaugh [l996[). 



7.4 Morphology 

In our model, bright elliptical galaxies form predominantly 
through galaxy mergers. When two galaxies of comparable 
mass coalesce (M2 > /ciupA/i, where Mi > M2), a violent 
merger is assumed to occur, leaving an elliptical galaxy as 
the remnant. Spheroids can also be built-up by the repeated 
accretion of smaller, gas-poor galaxies, because accreted 
stars are assumed to add to the bulge component of the ac- 
creting galaxy. Thus, the parameters fdf and /diip, which de- 
termine, respectively, the frequency of galaxy-galaxy merg- 
ers and the threshold above which a merger is deemed to 
be violent, are the primary parameters that influence the 
production of elliptical galaxies. The morphological mix de- 
pends also on the strength of stellar feedback. If feedback is 
weak, massive disks form at high redshift and have a long 
interval of time during which they can merge to form el- 
lipticals. Conversely, if feedback is strong, the formation 
of massive stellar disks is delayed, they experience fewer 
mergers and fewer ellipticals are produced. We can there- 
fore constrain these parameters by comparing the relative 



Table 3 . The morphological mix of galaxies brighter than Mb — 
51og/i < —19.5, for various values of the merger parameters /^f 
and /cllip and the feedback parameter Vhot- Also listed are two 
variants, which are described in the text. In the first, t unstable 
disks are transformed to spheroids and in the second,-!- accreted 
stars are added to the disk rather than the bulge. 



/cllip 


fdC 


Vhot/kms 1 


S : SO : E 


0.3 


1 


200 


61:08:31 


0.5 


1 


200 


70:07:23 


0.3 


0.5 


200 


53:09:38 


0.3 


2.0 


200 


79:06:15 


0.3 


1 


100 


38:05:57 


to.3 


1 


200 


46:09:45 


tO.3 


1 


200 


66:08:26 



abundances of galaxies of different morphological types in 
the model with observations. 

Our models do not strictly predict galaxy morphology, 
but rather the relative masses and luminosities of the bulge 
and disk components. The bulge-to-disk luminosity ratio is 
known to correlate with morphology, albeit with quite a 
large scatter, and so we simply take cuts in this ratio in 
order to define morphological classes. Ellipticals are defined 
as galaxies for which the bulge contributes more than 60% 
of the B-band light, spirals as those whose bulge contributes 
less than 40% of the B-band light, and SOs as galaxies in 
the intermediate range. The B-band magnitudes used here 
include dust extinction for galaxies with random inclination 
angles. 

Table 3 gives the predicted morphological mix of galax- 
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ies at the present day that resuhs from these definitions, for 
various values of the parameters fdf, /ciiip and Vhot- These 
ratios apply to a volume- limited sample of galaxies with ab- 
solute magnitude brighter than Mb = —19.5 + 51og/i, but 
the mix is very similar if one instead constructs an apparent 
magnitude-limited catalogue. For comparison, the morpho- 
logical mix in the APM Bright Galaxy Catalogue (which is 
apparent magn itude limited) is S -I- Irr : SO : E = 67 : 20 : 13 
(Loveday 1996, table 10), when one groups together spirals 



and irregulars, and assumes that the 90% of the galaxies 
in this survey that were classified are representative. This 
agrees well with the estimate S -I- Irr : SO -I- E = 66 : 34 for 
galaxies brighter than Mb = —19.0 + 51og/i i n the SSRS2 



redshift survey (Table 2 of Marzke et al. 1998 ). Comparing 



the model results with the observational values indicates 
that the morphological mix depends somewhat on the val- 
ues of Viiot , /ciiip and /df . Reasonable agreement with the ob- 
served mix can be obtained for a range of values of these pa- 
rameters, indicated by the examples given in the first three 
rows of Table 3. For the reference model (cf. Table 1), the 
values of /oiup and /df are approximately the lowest that 
would seem reasonable, given the physical processes that 
these parameters are trying to describe. Considering the 
crude manner in which we have defined morphologies the 
agreement between model and data is satisfactory, and it 
does not seem warranted to fine-tune these parameter val- 
ues any further. 

The morphological mix may al so be infiuenced by the 
disk instability discussed in Section 4.3.3. The penultimate 



row in Table 3 gives the mix found for the model with disk 
instability whose disk scalelength distribution was discussed 
in Section Here, it has been assumed that gas and stars 
in unstable disks are transferred to the bulge component, 
with the gas being consumed in a burst. We see that this 
significantly reduces the spiral fraction and boosts the el- 
liptical fraction, though the majority of ellipticals are still 
formed by mergers. In fact, the effect on the morphological 
mix is quite a strong function of luminosity. Brighter than 
Mb — —20.5 + 5 log h disk instability has very little effect, 
but it becomes increasingly important in low luminosity sys- 
tems. We plan to discuss the effects of disk instability in 
detail in a subsequent paper. 

One further assumption of our reference model is that 
stars that are accreted in minor mergers a re ad ded to the 
bulge of the resulting galaxies (see Secion 4.3.2). The last 
row in Table 3 shows how the morphological mix is modified 
if instead these accreted stars are added to the galaxy disks. 
Such a change only causes a modest increase in the spiral 
fraction. Also, we can see in Table 2 that this model (labelled 
"Accretion by disk") differs little from the reference model. 



7.5 Cold gas in Spiral Galaxies 

In our model, there are two parameters, e* and a*, which 
determin e the star formation timescale in galaxy disks (see 
equation 4.14). The first, e*, determines the star formation 



timescale for spiral galaxies with circular velocities compara- 
ble to those of Li, galaxies, while the second, a*, determines 
how this timescale varies with circular velocity. The lumi- 
nosity function (after rescaling by T) is fairly insensitive to 
e* and a*, but they do strongly affect the cold gas content 
of galaxies. 
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Figure 9. The cold gas content of spiral and irregular galax- 
ies as a function of luminosity. In each panel, the filled squares 
and associated errorbars show observational estimates of the me- 
dian, 10 and 90 percentile points respectively of the distribution 
of Miiy/Lg for Sa to Im galaxies. Here, Mjiy = Mm + is 
the mass of cold hydrogen which includes gas both in atomic and 
molecular forms. We have made these estimates using a coi nbina - 
tion of data from Huchtmeier & Richter ( |l988| ) and Sage ( |l993| ), 
as described in the text. The model results are shown by the 
open circles and their errorbars. For these, we select galaxies of 
comparable morphological type by requiring the B-band bulge-to- 
total luminosity ratios to be less than 0.4 . We express the model 
cold gas mass, M^^i^, in the observational units, h~^MQ, and 
set Mjjy = O.TMeoid to take account of the mass fraction of He. 
The upper panel shows the model with «^ = 0. The middle panel 
shows the reference model, which has Oi, = —1.5. The lower panel 
shows two models (the errorbars have been removed for clarity), 
each with a* = —1.5, but with the parameter e*, which controls 
the star formation timescale, varied up and down from the value 
in the reference model. 
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Fig. ^ shows how the amount of cold gas present in spi- 
ral and irregular galaxies depends on gala xy lum inosity. The 
observational data are taken from S age (|l993 ) and Hucht- 



meier &i Richter (1985). The Sage ( 1993[ ) data come from 
a complete sample of Sa-Sd galaxies with measurements of 
atomic HI and molecular H2, whereas the Huchtmeier & 



Richter (1988) data come from a complete sample of Sa- 



Im galaxies, but with only HI measurements. Brighter than 
AIb — 51og/i < —16, Sa-Sd galaxie s dom inate over Sdm-Im, 
and so we simply plot the Sage (199S) data. Fainter than 
Mb — 51og/i > —16, the mass fraction of molecular hydrogen 
appears to be small (Mhj/A^hi ^ 0.2), and so here we ne- 
glect the molecular H2 co ntribu tion and simply plot the data 
of Huchtmeier & Richter ( 1985 ). In each case, the luminosity 
is corrected to face-on. 

The model plotted in Fig. ^a has q* = 0, correspond- 
ing to the standard Kennicutt law in which the star forma- 
tion timescale is proportional to the disk dynamical time. 
In this case, the faint galaxies in the model typically con- 
tain less cold gas than is observed. In fact, the trend of 
Mgas/ie with Lb is in the opposite sense to that observed. 
The reference model plotted in Fig. |^, has q* = —1.5, 
implying longer star formation timescales for low circular 
velocity galaxies compared to a* = 0. This has a beneficial 
effect, and produces a correlation in Afgas/I/B vs. Lb which 
is much closer to the data. Fig. |^ shows the efi^ect of vary- 
ing e-^, and demonstrates how the observed gas fraction in 
bright spirals constrains this model parameter. 



7.6 Galaxy Metallicities 

Our model predictions for galaxy metallicities all scale lin- 
early with the yield p, aside from the effects of metallicity 
on the cooling of halo gas. We fix the value of p in our ref- 
erence model by requiring a good match to the mean stellar 
metallicity in L* ellipticals. 

In Fig. |l^, we show what the reference model predicts 
for the metallicity-luminosity relation for spiral and elliptical 
galaxies, compared to observational data. Fig. ^|a shows the 
gas metallicity vs. luminosity for spirals. The model points 
are derived from the metallicity of the cold, star-forming 
gas. T he observational data in this case, from Zaritsky et al. 
(1994), are based on HII region gas metallicities measured 
at r — 0.4_R25 in each galaxy (where R25 is the isophotal 
radius at a B-band surface brightness of 25 mag/arcsec^), 
rather than being averages over the whole galaxy. Fig. hOp 
shows the stellar metallicity vs. luminosity for ellipticals. 
The model points are obtained from mass-weighted mean 
stellar metallicities, but using V-band luminosity-weighted 
stellar metallicities would give nearly identical results. The 
observational data for the ellipticals are based on a compi- 
lation of estimates from line-strengths, which Zaritsky et al. 
have tried to put on a common metallicity scale. For both 
the spirals and ellipticals, we have converted the Zaritsky et 
al. observational data from relative to absolute metallicities 
assuming a solar metallicity Zq = 0.02. 

For both the spirals and ellipticals, our models predict a 
metallicity-luminosity relation which is in the same sense as 
the observed one, but is not as steep. The origin of this corre- 
lation lies in our model of stellar feedback. Our treatment of 
chemical evolution differs from the traditional "closed-box" 
model because gas reheated by SNe is allowed to escape from 
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Figure 10. The dependence of metallicity on luminosity in our 
reference model, compared with observational data. In each panel, 
the lines show the median metallicity in the model and the error 
bars indicate the 10 and 90 percentiles of the distribution. The 
obs ervati onal data, taken from the compilation by Zaritsky et 
al. (L994), are indicated by filled squares, where again the error 



bars show the 10 and 90 percentiles. The upper panel compares 
the metallicity of the cold star-forming gas in disk-dominated 
galaxies which in our model are galaxies selected to have B-band 
bulge-to-total ratios of less than 0.4, with observational data for 
spiral and irregular galaxies. The lower panel compares the stellar 
metallicity for bulge-dominated galaxies (B-band bulge-to-total 
ratios greater than 0.6) with observations for elliptical galaxies. 



the galaxy while cooling gas can fiow into it. Consequently, 
the appropriate expression for the effective yield becomes 
Pcff = (1 — e)p/(l — J?-|-/3) (equation B£), which depends on 
the s treng th of feedback via the quantity f3 given in equa- 
tion (4.15). Since f3 is smaller for larger, more massive galax- 



ies with deeper potential wells, their effective yield is large 
and this naturally results in more metal rich stellar popula- 
tions in more luminous galaxies. We could obtain a steeper 
metallicity-luminosity relation, in better agreement with the 
observed one, by assuming stronger feedback, but this would 
have deleterious effects on our fits to other properties. This 
issue is discussed further in Section 8.3 in connection with 
the colour-magnitude relation. 
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7.7 Mass-to-Light Ratios 

While the luminosities of galaxies are easily measured, the 
masses of the stellar populations they contain are less ac- 
curately determined, mainly because of the difficulties in 
separating the contributions to the mass from stars and 
dark matter in dynamical measurements. The mass-to-light 
ratios, M/L, for stellar populations in galaxies are corre- 
spondingly somewhat uncertain. For this reason, we do not 
use them as primary constraints in determining the model 
parameters. Nonetheless, they provide useful a consistency 
check, and can be used to exclude, for example, models 
which have very large brown dwarf fractions (i.e. large T), 
which might otherwise be viable. 

Table 2 lists the mass-to-light ratios of the stellar pop- 
ulations of both spiral and elliptical galaxies for each of our 
models. These mass-to-light ratios, which include the contri- 
bution of brown dwarfs, depend on the age and metallicity 
of the stellar populations, and also on the value of the pa- 
rameter T which has been fixed by reference t o th e bj-band 
luminosity function, as described in Section 7.1, Mass-to- 



light ratios are observed to depend on galaxy luminosity, 
so to make a fair comparison, we compare M/L values for 
observed and model galaxies at the same luminosity. The 
model M/L values given in the table are median values for 
—20 < Mb — 5 log h < 19. For each of the observational esti- 
mates described below, we have estimated the trend of M/L 
with luminosity from the values given for the individual 
galaxies in the paper concerned, and used this to estimate 
the average M/L for the galaxies at A/b — 51og/i ~ —19.5. 
For spiral galaxies, the observed values are dust-corrected 
to face-on values, and so the model M/L values are also 
calculated from face-on luminosities including the effects of 
dust. 

For spiral galaxies, most estimates of M/L are based 
on fit ting models to measured rotation curves. Buchhorn 
( [1992|) finds M/L = iAhMQ/L^ in the I-band, and Broeils 
(1992) 4.9/iM0/L0 in the B-band. These values are con- 
sistent with the mean colour, B-I ~ 1.8, found for spirals 



by de Jong ( 1996 ). Note, however, that these numbers are 
based on maximum disk fits to the observed rotation curves, 
and since a fraction of the rotation velocity may be pro- 
duced by dark matter, they should be viewed as upper lim- 
its to the stellar mass-to-light ratios. (For L* spiral galaxies 
in our reference model the mean contribution to the mass 
within the disk half-mass radius by non-baryonic dark mat- 
ter is 62%.) Comparing with the model values, we see that 
only the model with the high value of the baryon density 
(ilb = 0.04) comes close to violating these constraints. An 
alternative method for measuring the Ad/L of disks, which 
avoids contamination by dark matter, is to combine mea- 
surements of the vertical scaleheights and velocity disper- 
sions of galaxy disks. Using this method, Bottema (1997) 
finds an average M/L — 2AhMQ/LQ in the B-band, which 
is about a factor of two lower than the maximum disk value. 
The median M/L of our reference model and most of the 
variants are only slightly below this estimate, and so are 
quite compatible with observations. 

The comparison for el liptic al galaxies is similar. In th e 
B-band, Mobasher et al. ( [l999| ) and van der Marel ( |l99l| ) 
find M/L = 9.6hMQ /Lq and SHMq / Lq respectively, using 
stellar velocity dispersion measurements. Again, these values 



should be viewed as upper limits on the stellar mass-to-light 
ratios, as they include the effect of any dark matter within 
the effective radii of the galaxies. With the exception of the 
high baryon density model, all the models listed in Table 2 
are consistent with these data. 



7.8 Average Colours 

Table 2 also lists the median B-K colours of galaxies with 
—24.5 < Mk — 51og/i < —23.5 for various parameter val- 
ues. The observed median colour for this l uminosity r ange, 
calculated from the data of Gardner et al. (1996,1997) that 
are presented in Section is B-K — 3.8. The largest 
effects on the model colours result from varying the yield 
p, fib, and Viiot- A higher yield leads both to intrinsi- 
cally redder stellar populations and to greater amounts of 
dust, which redden the observed galaxy colours still further. 
Increasing fl^, increases the gas density in halos, and thus 
shortens the cooling time. This then results in a greater frac- 
tion of the gas cooling and forming stars at high redshift, 
producing an older, redder stellar population. Decreasing 
Vhot reduces the strength of feedback, allowing more stars 
to form in low-mass halos at high redshift, leading, again, to 
older, redder stellar populations by the present. Somewhat 
surprisin gly, i ncreasing the star formation efficiency, e*, (cf. 



equation 4.14) has little effect on the present star formation 



rate and galaxy colours. This happens because the shorter 
star formation timescale allows more star formation to occur 
at high redshift, and this leads to a reduction in the amount 
of cold star-forming gas at low redshift which compensates 
for the shortened star formation timescale. 

The colours also depend on the choice of IMF. However, 
changing from the Kennicutt to the Salpeter form makes 
very little difference to the B-K colours, which become bluer 
by just 0.015 magnitudes. This is to be expected, as Fig. ^ 
demonstrates that the differences in the spectral properties 
of the stellar populations for these two IMFs are small. The 
metallicities of the stellar populations and the dust con- 
tent of the galaxies are very similar in the two cases, as the 
adopted values of the yield are identical and of the recycled 
fraction almost identical. 



7.9 Summary of Parameter Constraints 

In this section we have demonstrated how the predictions 
of local galaxy properties depend on each of the model pa- 
rameters. We now summarise the reasons for adopting the 
specific parameter values that define our standard or refer- 
ence ACDM model. 

The cosmological parameters, Qo ~ 0.3 and Aq = 0.7, 
were chosen without reference to galaxy properties, and sim- 
ply define the background cosmology in which we sought a 
viable model of galaxy formation. The Hubble parameter, 
h = 0.7, was chosen to be in reasonable accord with recent 
estimates. The baryon density, fib, was chosen as a compro- 
mise between constraints from primordial nucleosynthesis 
and the need to prevent the disks of bright galaxies becom- 
ing strongly self-gravitating and so distorting the bright end 
of the TuUy-Fisher relation. This choice also prevents galax- 
ies becoming too luminous, or equivalently, their M/L ratios 
becoming too large once T has been adjusted. The shape. 
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r, of the mass fluctuation power spectrum was chosen to 
be c ons istent with the above choices of fio, and h (cf. 
eqn 7.1). The resulting value, F = 0.19, is in accord with 



the shape of the power spectrum inferred from studies of 
large-scale galaxy clustering. The amplitude, erg, was cho- 
sen for c onsistency with t he observed abundance of galaxy 
clusters (Eke et al. 1996). For our chosen cosmology, this 



is consistent with the value of as preferred by the COBE 
measurements of microwave background anisotropics. 

The stellar population parameters, we have essentially 
chosen to be consistent with models constructed to account 
for solar neighbourhood data. In particular, we have adopted 
the Kennicutt IMF. The fraction of gas recycled via stellar 
winds and SNe, we have taken to be i? = 0.42/T — 0.31, 
consistent with what is expected for this IMF. (Recall that 
T is defined as the total mass in stars, including brown 
dwarfs, divided by the mass in luminous stars. The value 
T = 1.38 was chosen by fitting the position of the bright 
end of the bj-band luminosity function.) We have adopted a 
yield, p = 0.02, which results in roughly the observed metal- 
licity for galaxies similar to the Milky Way and for L-^ ellipti- 
cals. This may be seen in Fig. ^ where we also examine the 
model prediction for the dependence of metallicity on galaxy 
luminosity. Our adopted yield implies pi = pT = 0.028, 
which is also roughly consistent with theoretical estimates 
for our assumed IMF. The yield also determines the metallic- 
ity and dust content of the cold gas disk. Our adopted value 
gives rise to an amount of reddening which is approximately 
that required to bring the model into good agreement with 
the observed galaxy B~K colour distribution. Varying the 
IMF between the Kennicutt, Salpeter or Miller-Scalo forms 
has only a small effect on the galaxy colours and mass-to- 
light ratios at 2: = 0, and so the IMF is therefore not well 
constrained by the data we have examined so far. However, 
the small differences in these IMFs can have a significant 
effect on model predictions at high redshift. 

The dynamical friction parameter, fdf, we simply set to 
its natural value, fdf = 1. 

The parameter fdup, which sets the threshold for vio- 
lent gal axy m ergers which produce ellipticals and bulges (cf. 
Section 4.3.2| ), was chosen so as to reproduce an acceptable 
mix of morphological types. 

Finally, the model requires parameters for the star for- 



mation and feedback laws (cf § 1.2). The feedback parame- 
ters, Vhot and cthot , we have constrained by the shape of the 
bj-band luminosity function, the slope of the Tully-Fisher 
relation, and the distribution of spiral disk scalelengths. A 
low value of Ohot is required to avoid curvature in the Tully- 
Fisher relation, while a large value helps to reduce the faint 
end slope of the galaxy luminosity function. Our adopted 
compromise, ai^ot ~ 2, results in a straight Tully-Fisher rela- 
tion and a luminosity function with a faint end slope in good 
agreement with the ESP survey. This is steeper than insev- 



eral other surveys, including those by Loveday et al. (1992) 



and Ratcliffe et al. (1998), but we have demonstrated that 



the slope is sensitive to surface brightness selection limits. 
The value of Vhot influences both the faint end of the lumi- 
nosity function and the sizes of galaxies. Our adopted value 
of Viiot = 200kms~^ appears to be a good compromise. The 
parameter, e, which allows metals produced in SNe to es- 
cape directly to the surrounding hot halo gas rather than 
flrst being mixed with the cold gas in the galaxy disk, we 



have simpl y se t to zero . Regarding the star formation law 
(equations 4.4 and 4.14), the overall scaling of the timescale, 
set by the efficiency e*, is constrained primarily by the cold 
gas masses in spirals. The dependence of this timescale on 
galaxy circular velocity is determined by a*. We adopted the 
value Of* = —1.5 in order to improve the correspondence be- 
tween model and observed cold gas masses in low-luminosity 
disk galaxies. 



8 FURTHER PROPERTIES OF THE 
REFERENCE MODEL 

In this section we compare further properties of our refer- 
ence model with observational data that have not been used 
to set the model parameters. We defer comparisons with ob- 
servations at high and moderate redshift to future papers. 

8.1 Colour Distributions 

We have already consider ed t he average B-K colours of 
L ^ Li, galaxies in Section 



7.S. In our reference model, the 



mean galaxy B-K colour is quite strongly constrained by 
virtue of the requirement that the model give a reasonable 
flt to the bright end of both the bj and K-band luminos- 
ity functions. Nevertheless, it is still interesting to look at 
the full distribution of predicted colours, and to compare 
them to other observational data. In Fig. |l^, the reference 
model is compared to the distributions of B-K and B-V 
colours in the K-selected redshift survey of Gardner et al. 
( [I996| , [i997| ). This contains more than 500 galaxies, covers 
10 square degrees and was imaged in the B, V, I and K 
bands. The redshift and colour information allow accurate 
k-corrections to be derived by matching each galaxy's ob- 
served colours with one of a set of template spectra. Thus, 
the histograms in Fig. ^ show rest-frame colours. If the 
more uncertain correction for evolution is also applied, then 
the observed distributions typically shift redwards by 0.15 
magnitudes. The model with dust matches both the median 
colours and the widths of the colour distributions quite well. 
If the effects of dust are not taken into account, then the re- 
sulting model median colours are too blue by approximately 
0.3 magnitudes in B-K and approximately 0.1 magnitudes 
in B-V . Thus, we see again how modelling the effects of 
dust is an important factor in producing acceptable galaxy 
colours. 



8.2 The Colour-Morphology Relation 

Fig. |l2| shows the predicted correlation of galaxy B-V 
colour, for face-on inclination, with bulge-to-disk ratio. A 
strong correlation is predicted, with bulge-dominated galax- 
ies typically being much redder than disk-dominated galax- 
ies. However, galaxies which have experienced a major 
merger and the associated burst of star formation within 
the last 1.0 Gyr have large bulge-to-disk ratios, but are 
much bluer than bulge-dominated galaxies that have not had 
a recent major merger. Observationally, "normal" galaxies 
show a similar trend in colour to the model galaxies which 
have not had recent bursts, as is illustrated by the obser- 
vational data of Buta et al. (1994) which are also plotted, 
with the solid line showing the mean galaxy colour and the 
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Figure 11. A comparison of galaxy colours in the reference model with observations. The upper two panels compare rest-frame B-K 
colour distributions in volume-limited samples for two different ranges of absolute K-hand magnitude. The lower two panels compare the 
rest-frame B—V colour distributions for the same two ranges of absolute K-hand magnitude. In e ach case, th e histograms with errorbars 
show the observational distributions, derived from the K-haad redshift survey of Gardner et al. (1996 1997) using the l/Vmax method. 
These data have been k-corrected by matching each galaxy's observed colours with one of a set of template spectra. If the more uncertain 
correction for evolution is also applied, then the observed distributions typically shift redwards by 0.15 magnitudes. The dotted lines show 
the colour distribution of the reference model without including the effects of dust and the solid lines show the distributions including 
the effects of dust. 



dotted lines showing the dispersion. Note that Buta et aL 
have removed from their sample some galaxies viewed as 
outliers or having strong emission lines. To plot these data 
on Fig. we have converted the Hubble T-type given in 
Table 6 of Buta et al. to bulge-to-total ratio, using the fi t 
given in equation (5) of Simien & de Vaucouleurs (1986). 
This fit represents an extrapolation from T = —3 (which 
corresponds to B /T ~ 0.6) to T = —5 (which corresponds 
to B/T « 1.0). Given the considerable scatter that exists 



between bulge- to-tota l ratio and T-type (e.g. see Figure 1. 
of Baugh et al. |l996b| ), the level of agreement between the 
predictions of the model and the data is encouraging. It will 
be interesting to perform more quantitative comparisons of 
this prediction with observations when a suitable data set 
exists in which bulge-to-disk decompositions have been done 
for a complete survey. 
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Figure 12. Average galaxy B-V colour as a function of bulge- 
to-disk ratio, compared to observations. The open squares and 
error bars show the model mean B-V colour and rms dispersion 
as a function of B-band bulge-to-total ratio for galaxies brighter 
than My — 5 log ft. = —20.5. Face-on colours are plotted, including 
the effects of extinction. The solid square shows the mean colour 
of model galaxies which have experienced a merger and the as- 
sociated burst of star formation in the last 1 Gyr. The solid and 
dotted lines show the observed mean and rms dispersion of the 
colour from the data of Buta et al. (1994). Hubble T-type has 
been converted to bulge-to-total ratio using the data of Simien & 
de Vaucouleurs (1986), as described in Baugh et al. (1996b). 
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Figure 13. The colour-magnitude relation for cluster elliptical 
galaxies in the reference model, compared to observations. The 
points give the predicted distribution of V-K colour versus V-band 
magnitude for elliptical galaxies in clusters with circular velocity 
greater than 1000km s~^. The heavy line and error bars indicate 
the median and the 20 and 80 percentiles of this distribution. 
The o bserved correlation and scatter, from Bower, Lucey & Ellis 
(1995), are indicated by the dotted line and associated error bars. 



to a flattening in the predicted colour-magnitude relation 
for bright ellipticals. Our models are capable of producing a 
significant gradient in this diagram if we assume very strong 
f eedbac k and a large yield p, just as Kaufi'mann & Chariot 
(1998a) did, but this is then at the expense of other suc- 
cesses of the model. In particular, strong feedback gives rise 
to excessively large disk scalelengths. This is clearly an issue 
that deserves further investigation. 



8.3 The Elliptical Colour-Magnitude Relation 

A second correlation that we have examined is the colour- 
magnitude relation of elliptical galaxies. This is closely 
related to the metallicity-luminosity re latio n for elliptical 
galaxies, 



7.6 



Fig. |l| 



aheady discussed in Section 
pares the distribution that we predict for cluster ellipticals 
with that determined o bserv ationally for the Coma cluster 
(Bower, Lucey & Ellis 1992). The model predicts quite a 
small spread in the colours of bright ellipticals, consistent 
with the observed scatter, but it does not reproduce the 
strength of the observed correlation of colour with luminos- 
ity. This is the same problem that we fou nd when we made 
this comparison in Baugh et al. (1996b), even though our 



new model includes chemical enrichment, which the earlier 
model did not. 



In Kau ffman n et al. (1993), Baugh et al. (1996b) and 
Kauffmann ( ll99^ ), it was argued that the inclusion of chem- 
ical enrichment (which was neglected in the early models) 
might give rise to the desired correlation. This was explic- 
itly demonstrated for certain models by Kauffmann & Char- 
lot (1998a). In our present models, chemical enrichment 



does appear to have the desired effect at low luminosities, 
where the models produce a similar gradient in the colour- 
magnitude relation to the one observed. However, the cor- 
relation between metallicity and luminosity flattens at the 
brightest magnitudes (see Figure |lOp), and this gives rise 



8.4 The cosmic star formation history 

It is interesting to examine how the improvements in our 
modelling techniques affect our predictions for the globa l 
history of star formation, first presented in Cole et al. ( 1994 ) . 
Fig. hj shows the redshift evolution of the mass fractions of 
baryons in the forms of hot gas, cold gas and stars, and the 
mean metallicities of each of these components. Consistent 
with our feedback model, we have assumed that baryons con- 
tained in halos with mass below our resolution limit are in 
the hot, diffuse phase. We remind the reader that by "hot" 
gas we simply mean diffuse gas in halos, whatever its physi- 
cal temperature, and by "cold" gas we mean all the gas that 
has cooled and collapsed into galaxies. The mass of stars 
increases steadily towards the present, with just over half of 
the present stellar mass having been formed since redshift 
z < 1.5. This is similar to our earlier results in Cole et al. 



(1994) and Baugh et al. (1998), despite the various changes 
in model parameters that we discuss below. The evolution 
of the cold gas mass shows a broad peak between redshifts 
1 < z < 2. The fall-off at high redshift reflects the efficiency 
of stellar feedback, which keeps gas in low mass halos in the 
hot phase. At low redshift, the cold gas fraction begins to 
decline as cooling becomes increasingly less efficient in the 
large mass, high virial temperature groups and clusters that 
form at late times, while, at the same time, the reservoirs of 
cold gas are depleted by ongoing star formation. The mean 
metallicity of the hot gas grows at a rate which mirrors the 
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Figure 15. The luminous star formation rate (i.e. excluding 
brown dwarfs) per unit comoving volume as a function of redshift. 
The solid line shows the total star formation rate per unit volume 
in the reference model. The dashed line shows the contribution 
from quiescent star formation in galactic disks. The dotted line 
shows the contribution from bursts of star formation that occur 
during major mergers. The dot-dashed line is the star formation 
history in model G of Baugh et al. (1998) (for the same cosmologi- 
cal parameters as in our reference model) . In the current reference 
model, the star formation rate per unit volume is higher at high 
redshift than in the old model, because of the weaker feedback in 
low mass halos and the shorter star formation timescale at high 
redshift that are assumed in the new model. 



Figure 14. Evolution of baryonic mass fractions and average 
metallicities. The upper panel shows the fraction of baryons in 
stars, cold gas and hot g function of redshift, and the 

lower panel shows their corresponding mean metallicities. The 
solid lines are for stars, the dashed lines for cold gas and the 
dotted lines for hot gas. 



growth in total stellar mass. In contrast, the mean metallic- 
ity of the stars and cold gas reaches 1/3 of its present value 
even at very high redshift, and then increases only gradually 
between redshift 2 = 6 and 2 = 0. 

Fig. |l^ shows the star formation history in differential 
form. The solid line is the average formation rate of luminous 
stars per unit comoving volume (i.e. the total star formation 
rate divided by T). The dashed and dotted lines show sepa- 
rately the contributions to the total rate from quiescent star 
formation in disks and from bursts of star formation induced 
by galaxy mergers. Approximately 10% of the stars formed 
at any redshift are formed in bursts. 

The dot-dashed line in Fig. |l^ shows our earlier result, 
presented as model G of Baugh et al. (1998). Model G had 



very similar cosmological parameters to the present refer- 
ence model, but was calculated with a slightly earlier ver- 
sion of our code. The differences between the two results 
are easy to understand and provide a nice illustration of 
how various aspects of our galaxy formation modelling are 
closely interconnected. The main difference between the two 
model predictions is the lower star formation rate at high 
redshift obtained in model G compared to the present model. 
This difference mainly reflects the different values of the 
feedback parameters that we adopted in the two models. 
In our older model we set these parameters by the require- 
ment that the faint-end of the present-day galaxy luminosity 



function reprod uce a s closely as possible that measured by 
Loveday et al. 01992 ). This led us in Baugh et al. (and in 
Cole et al. 1994|) to assume very strong feedback. Here, we 



have argued that the uncertainty in the faint end of the 
luminosity function as indicated by the wide range of ob- 
servational esima tes, i ncluding the very steep slope found 
by Zucca et al. (1997), suggests that this is not a robust 
constraint. We have instead adopted a weaker feedback that 
avoids introducing curvature in the faint end of the TuUy- 
Fisher relation. Thus, in our models, the behaviour of the 
star formation rate at high redshift is intimately tied in with 
our assumptions about the strength of stellar feedback and 
the faint-end of the present-day galaxy luminosity function. 
A second diffe rence between the present model and that of 
Baugh et al. ( [1998| ) is the assumed dependence of the star 
formation timescale on galaxy properties. We now incorpo- 
rate a scaling with the galaxy dynamical time, which results 
in all star formation timescales at high redshift being smaller 
than in the older model. We must emphasize, however, that 
in spite of the uncertainties in the predicted star formation 
rate at high redshift, our prediction of a late epoch for the 
m ajority of star formation is robust ( cf. Fig ure 21 of Cole et 
al 



1994 



and Figure 14 of Baugh et al. 1998| ) . This is because 
in all the versions of our model, only a small fraction of stars 
form at 2 > 3, and it remains the case that we expect half 
the stars in the universe to have formed at redshift 2 ^ 1.5. 



9 DISCUSSION 

In this paper, we have presented a new semi-analytic model 
of gal axy fo rmation, based upon the one developed by Cole 
et al. (1994). Our new model contains a number of additions 
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and improvements. For example, we have designed and im- 
plemented a new algorithm to generate halo merger trees 
with arbitrary mass resolution, and we have extended our 
modelling to include realistic descriptions of the density pro- 
files of dark matter halos and their gas content, as well as cal- 
culations of chemical evolution, dust extinction, and galaxy 
sizes. We applied this model to a specific cosmology, the 
ACDM model, which has fio = 0.3, Ao — 0.7 and a pri- 
mordial power spectrum whose amplitude is consistent with 
both the local abundance of galaxy clusters and with the 
COBE anisotropics in the microwave background radiation. 
We then compared the results with a range of observational 
data for the local galaxy population. 

In principle, a model like ours can predict virtually any 
simple property of the galaxy population over a large range 
of redshift. However, not all predictions are equally reliable. 
For a given cosmological model (e.g. ACDM), the evolu- 
tion of the population of dark matter halos is known with 
high accuracy and can be calculated without any free pa- 
rameters. The initial internal structure of these halos is also 
completely specified if one adopt s the results of recent high- 
resolution N-body simulations (Navarro et al. 1997). The 
process of gas cooling is more uncertain but can be calcu- 
lated also without free parameters (other than the value of 
the gas metallicity), using a model based on simple assump- 
tions about the geometry and initial configuration of the gas. 
We assume an initial spherically symmetric distribution of 
gas at the virial temperature of the halo in which it is con- 
tained, and a gas density profile given by the "/3-model". 
In practice, the dynamics of the cooling gas are likely to 
be subtantially more complex than this s imple model im- 
plies. Nevertheless, as Benson et al. ( 2000c ) have shown, the 



simple model does predict global fractions of hot and cold 
gas, and their distribution in halos of different mass, that 
are broadly in agreement with the results of gasdynamics 
simulations. 

The formation of stars from the gas that has cooled 
and the associated feedback processes are the m ost un cer- 
tain components of the model. As in Cole et al. (1994), we 
have represented them using simple scaling laws, but we 
have adopted a more flexible treatment than we had done 
previously. In addition to specifying the IMF and the asso- 
ciated fraction of brown dwarfs, our model of star formation 
and feedback requires four free parameters. Our treatment 
of feedback is simplified and neglects potentially important 
sources of energy such as active galactic nuclei and quasars. 
It also neglects the dynamical response of the hot halo gas 
to the heating generated by the injection of feedback en- 
ergy. In principle, this energy can modify the structure of 
the gaseous halo and hence its cooling rate. Thus, the de- 
tailed treatment of feedback impacts on all aspects of the 
galaxy formation model. For example, as we discussed in 
Section 7, our model only works well if we assume a value 
of the mean cosmic baryon density, Qi,, which is lower than 
recent determinations based on the deuterium abundance 



at high redsh ift by Schramm & Turner (1998) and Buries & 
Tytler (199J), although it is co nsistent with older dete rmi- 



nations by Walker et al. (1991) and Copi et al. (1995). A 



larger value of fib causes too much gas to cool, leading to a 
poor match to the TuUy-Fisher relation and to unacceptably 
large mass-to-light ratios. A successful model with larger val- 
ues of flb, however, is likely to be possible if feedback raises 



the entropy of the hot halo gas, thereby stron gly su ppress- 
ing cooling, as argued recently by Bower et al. (200C) in the 
context of X-ray clusters. 

Intimately linked to the processes of star formation and 
feedback is the chemical evolution of the gas. Although the 
basic principles of chemical enrichment are well-understood, 
important aspects, such as the mixing of metals in the in- 
terstellar medium, remain uncertain. Our model of chemical 
evolution requires specifying 3 parameters, two of which, 
however, (the yield and the fraction of stellar mass ejected 
by stellar winds and SNe) are determined by the choice of 
IMF. The remaining parameter is related to the mixing of 
metals. Once the metallicity of the gas is obtained, the spec- 
trophotometric properties of the stars are calculated using 
a population synthesis model which has no free parameters. 

These five ingredients: gravitational evolution of halos, 
gas cooling, star formation and feedback, chemical evolu- 
tion, and stellar population synthesis make up the core of 
our model of galaxy formation. To predict observable quan- 
tities, however, still requires a model for dust extinction. We 
assume that the mass of dust that forms is proportional to 
the mass in metals, and derive the extinction in any pass- 
band using a simple model for the distribution of dust in the 
disk. Our dust model has one free parameter (the ratio of 
the scaleheights of dust and stars), but our results are insen- 
sitive to its value. The core galaxy formation model can now 
be used to predict a wide range of visible galaxy properties 
including basic quantities such as the luminosity function in 
different passbands or the distribution of colours, as well as 
their evolution with redshift. 

To go beyond estimates of luminosity requires the ad- 
dition of more physical ingredients into the model. As the 
model becomes increasingly complex, it encompasses a fuller 
range of galaxy properties and this, inevitably, requires a 
growing number of physical assumptions and additional pa- 
rameters. For example, it is possible to distinguish between 
disk and spheroidal stellar conflgurations by introducing 
simple but plausible assumptions. Here we have assumed 
that cooling gas settles into a centrifugally supported disk, 
that the distribution of halo and gas angular momentum 
has a particular form (consistent with results of N-body 
simulations), and that major mergers or disk instabilities 
produce spheroidal stellar systems. Our model requires one 
further free parameter to deflne what a major merger is. 
With these assumptions, the scalelengths of disks can be 
computed without further free parameters, as can the sizes 
of spheroids, assuming that energy is conserved in mergers. 

In summary, a model of galaxy formation consists of 
a mixture of assumptions about the physical processes at 
work, together with adjustable parameters that reflect our 
lack of knowledge about certain complex astrophysical pro- 
cesses such as star formation. It is important to recognize 
that these parameters are not statistical variables describing 
a particular dataset (e.g. the faint-end slope of the luminos- 
ity function), but genuine physical quantities that describe 
a model for a speciflc physical process (e.g. the conversion of 
cold gas into stars). In many instances, the parameters can 
vary only over a relatively narrow range of physically sensi- 
ble values. In our model, just as in models by other groups, 
we strive to make the simplest possible assumptions at every 
stage and to introduce the minimum number of adjustable 
parameters, the majority of which, in fact, are required to 
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describe the poorly understood processes of star formation, 
feedback and the mixing of metals. Given the current theo- 
retical and observational understanding of these processes, it 
is not possible to build a realistic model of galaxy formation 
which has substantially fewer parameters than ours. 

Although the number of parameters is small, in any 
case, compared to the vast array of properties that the model 
can predict, it is important to adopt a well-defined, a priori 
methodology for fixing their values and for testing the va- 
lidity of physical assumptions. In our case, this strategy is 
straightforward: we fix the values of all the parameters by 
attempting to match a small subset of the local galaxy data 
which we regard as the most fundamental. In order of the 
weight we give them, these are: the luminosity functions in 
the B- and K-bands; the relative fractions of ellipticals, SOs, 
and spirals; the slope of the TuUy-Fisher relation at faint 
magnitudes; gas fractions in disks as a function of B-band 
luminosity; the distribution of disk sizes; and the metallicity 
of L* ellipticals. Once the model parameters have been set 
by these comparisons, we test our model predictions against 
a wide range of other data, without any further adjustments. 

The galaxy formation model presented here dif fers i n 
several significant respects from that of Cole et al. (1994). 
That model was based on a cruder method for generating 
halo merger trees, but that alone makes little difference to 
the resulting galaxy properties. Similarly, the extensions we 
have included which allow us to predict more galaxy proper- 
ties, such as sizes and mean metallicities, make little differ- 
ence to the properties that we were able previously to pre- 
dict. There are three differences which do lead to significant 
changes in our results. Firstly, the adoption of a cosmolog- 
ical model with a low value of f2o reduces the number of 
galaxy-sized halos, and this helps to reduce the offset in the 
TuUy-Fisher relation for models normali zed to the B-band 
luminosity function (see Heyl et al. 1995). Secondly, the in- 



clusion of dust in the calculation of luminosities and colours 
makes a typical galaxy colour slightly redder and this helps 
to match the B and K-band luminosity functions simulta- 
neously. Furthermore, since the luminosities that enter into 
the I-band TuUy-Fisher relation are partially corrected for 
the effects of extinction, the inclusion of dust also helps re- 
duce the offset between model and data seen in Cole et al. 



(1994). Thirdly, we have adopted a weaker feedback law for 
low circular velocity galaxies, since some recent determina- 
tions of the galaxy luminosity function indicate that the very 
flat faint-end which we strived hard to match in Cole et al. 



(1994) is not a robust observational constraint and may be 
affected by survey selection criteria. This, in turn, implies a 
higher star formation rate at redshifts z <: 3. 

Although in this paper we have focussed on the ACDM 
model, we have also investigated models with other values 
of the cosmological parameters. For reasons of space, we 
do not present our results in any detail here, but confine 
ourselves instead to a few general remarks, applicable to 
cluster-normalized models. A standard CDM model (SCDM, 
J7o = 1) still fails to match the TuUy-Fisher relation (even 
using halo circular velocities) and the luminosity function 
simultaneously. An f2o = 1 model with the same power 
spectrum shape as ACDM (the rCDM model of Jenkins et 
al. 1998) shares this problem and, in addition, the smaller 



both a lower star formation rate density and a lower abun- 
dance of Lyman-break galaxies at high redshift. (This prob- 
lem does not afilict the ACDM model because the effect due 
to the shape of the power spectrum is compensated for by 
the difference in the linear growth rate and the higher initial 
amplitude required by the cluster normalization.) 

There are now in the literature a number of semi- 
analytic galaxy formation models of varying sophistication, 
based on s imilar principles to ou rs (e.g. Avila-Rees & Fir- 
Guiderdoni et al. 



mani 1998 



1995, Roukema 1998, Wu et al 



199|, Kauffmann et al. |l999^SomerviUe & Primack |999|) 



It is beyond the scope of this paper to carry out a detailed 
comparison of all these models, but we re fer th e reader to 
the recent paper by Somerville & Primack (199£) which has 
compared some of the different approaches, including those 
of the Durham, Munich, and Santa-Cruz groups. For the 
most part, results from different models tend to agree well 
when similar assumptions are made. In practice, however, 
it is not uncommon for different groups to make somewhat 
different assumptions and, most importantly, to include dif- 
ferent physical effects in their models. This naturally leads to 
different results. An example of the former are the different 
strategies for constraining the star f ormatio n and feedback 
laws adopted by Kauffmann et al. (1999a) and ourselves. 
We give most weight to the local B-band luminosity func- 
tion and do not make any further adjustments when calcu- 
lating the zero-po int of TuUy-Fisher relation. By contrast, 
Kauffmann et al. (1999a) give most weight to the zero-point 
of the TuUy-Fisher relation. Since the models do not match 
both these observables perfectly, there is a difference in the 
luminosity normalization of the two models that propagates 
to other observables. An example of different physical pro- 
cesses is the treatment of the structure and angular momen- 
tum transport of gas cooling onto galaxies adopted by Wu et 
al. ( |1998D . Their model leads to a completely different mech- 
anism for making ellipticals and spirals to the one operat ing 
in our own model or in that of Kauffmann et al. (1999a). 



10 CONCLUSIONS 

We have presented a new semi-analytic model of galaxy for- 
mation which contains several novel features. It employs a 
state-of-the art Monte-Carlo algorithm for calculating the 
merging evolution of dark matter halos and it incorporates, 
for the first time, detailed prescriptions for calculating the 
sizes of disks and spheroids. We used this model to calculate 
observable properties of galaxies in the ACDM cosmology 
(fio = 0.3, Ao = 0.7, h = 0.7, and ug = 0.93) and focussed 
primarily on galaxy properties at the current epoch, with 
the following main conclusions: 

1. A pleasing agreement can now be obtained between the 
model and observed galaxy luminosity functions in the B- 
band and the K-band, over at least 8 magnitudes. This is 
a non-trivial success. In the B-band, the model was tuned 
to fit the ESP luminosity function of Zucca et al. (1997), 



amount of small-scale power compared to SCDM results in 
a later epoch for the onset of galaxy formation and, thus, in 



which has a steep faint-end. Unfortunately, there is still a 
large uncertainty in the observational estimate of the num- 
ber of galaxies fainter than L*. This is disappointing be- 
cause the faint-end slope of the luminosity function is ex- 
tremely sensitive to feedback processes which are therefore 
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only crudely constrained. A flatter fai nt-en d, like that mea- 
sured, for example, by Loveday et al. (1992) in the Stromlo- 



APM survey, could be obtained by increasing the strength 
of feedback. Our inability to constrain the feedback model 
better has a knock-on effect on our ability to predict the 
cosmic star formation rate at high redshift since this too is 
strongly influenced by feedback. Dust extinction has a rela- 
tively modest effect, dimming the bright end of the B-band 
luminosity function by about 0.5 mag. We showed that sur- 
face brightness effects can be important for faint galaxies, 
and this could help explain some of the discordant estimates 
of the faint end of the luminosity function. 

2. Our model reproduces both the observed mean galaxy 
colours and the spread in colour over a large range of galaxy 
luminosity. The stellar mass-to-light ratios of both stellar 
disks and ellipticals match the observational values well. In- 
clusion of reddening is important for this comparison, which 
does not involve adjusting any model parameters. The mean 
and scatter in the colours of galaxies of different morpholog- 
ical types, as measured by blue bulge-to-total light ratio, are 
also reproduced well. 

3. The current cold gas content of galaxies of different lu- 
minosity is related to the efficiency of past and current star 
formation. Our adopted star formation model (which iscon- 
sistent with the observations analysed by Kennicutt 199J) 



leads to excellent agreement with the observed ratio of cold 
gas mass to blue luminosity over 7 magnitudes. 

4. The predicted distribution of disk sizes is sensitive to 
the strength of feedback. Our model agrees well with the 
data, particularly for bright galaxies. 

5. The more realistic treatment of various properties and 
processes (e.g. dark halo and gas density profiles, dust, etc) 
leads to a better match to the I-band Tully-Fisher relation 
than was possible with our earlier, simpler model. If, as in 
most previous work of this kind, the circular velocity of a 
galaxy is identified with the circular velocity at the virial 
radius of the halo in which it formed, then our model gives 
an excellent fit to the zero-point, slope and scatter of the 
Tully-Fisher relation. However, in the model, the rotation 
velocities of galaxy disks at their half-mass radii are typ- 
ically 30% higher than the circular velocities of the halos 
at the virial radius. This results in an offset of -1-30% in the 
velocity zero-point of the Tully-Fisher relation when the cal- 
culated disk velocities are used. It remains unclear whether 
this disagreement reflects a fundamental shortcoming of the 
cold dark matter theory or whether it is simply a reflec- 
tion of various physical uncertainties in the calculation. For 
example, the derived disk rotation velocity depends on the 
assumptions of angular momentum conservation and adia- 
batic invariance during the collapse and formation of a galac- 
tic disk. If the collapse were, in fact, clumpy, then angular 
momentum would be transferred fro m th e disk to the halo 
(Frenket al. |l984 Navarro & White [f994 Navarro & Stein- 
metz 1997). In this case our calculation may have overesti- 



mated the amount by which the inner part of the halo con- 
tracts. This would help reduce the Tully-Fisher offset, but at 
the same time the loss of angular momentum from the disk 
would make the disks physically smaller and could act in 
the opposite direction, compressing the halo more strongly. 
These issues are worthy of further investigation, but they are 
best ad dressed usin g numerical simulations (see e.g. Navarre 
& Steiijmetz 1999|). 



6. Our model calculates chemical evolution, taking into ac- 
count the effects of gas loss due to winds and gas accretion 
due to cooling in a self-consistent way. The model predicts 
a trend of increasing metallicity with luminosity, similar to 
that observed, for star-forming gas in disk-dominated galax- 
ies and for stars in bulge-dominated galaxies. However, the 
colour-magnitude relation for ellipticals in clusters is signifi- 
cantly flatter than observed at bright magnitudes, al though 
the scatter is about right. Kauffmann & Chariot (1998e) 



have shown that a steeper slope for the colour-magnitude 
relation can be obtained by simultaneously increasing the 
strength of the feedback and the value of the yield. How- 
ever, in our ACDM model, a much stronger feedback than 
we have assumed would result in disk sizes that are much too 
large (because the accretion of gas onto galaxies is delayed), 
and would degrade the fit to the Tully-Fisher relation. We 
intend to carry out a more thorough investigation of this 
conflict in a later paper. 

7. Our more sophisticated modelling te chniqu es do not 
ch ange our earlier conclusion (Cole et al. 1994, Baugh et 
al. 1995) that half of the stars in the universe formed since 



2 ^ 1.5. However, the relaxation of the requirement for 
strong feedback (arising from the fact that we now fit the 
steep faint-end slope of the ESP luminosity function rather 
than the flat slope of the Stromlo-APM survey) allows a 
somewhat higher star formation rate a,t z ^ 3 than we had 
predicted previously. The fraction of baryons in cold gas has 
a broad peak at 1 < z < 2. The evolution of the mean metal- 
licity of the hot gas mirrors t he gr owth of stellar mass but, 
as noted also by Kauffmann (1996), the mean metallicity of 



the stars and cold gas builds up very rapidly: it is already 
about one third of the present value at z = 5. 

In summary, the model we have presented is broadly 
successful in matching a large range of galaxy properties. 
There remain, however, some interesting discrepancies, for 
example, the Tully-Fisher relation and the colour-magnitude 
relation for cluster ellipticals. Although the discrepancies are 
relatively small, further work is required to assess whether 
they point to incorrect assumptions or to the neglect of im- 
portant physical processes in our modelling procedure. 
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APPENDIX A: HALO ROTATION VELOCITY 

In this appendix we relate the halo rotation velocity, Kot, 
(assumed constant) to its spin parameter Ah. 

The total angular momentum of the halo is given by 



that r^p(r)(T^(r) vanishes as r — > we obtain, 



— Kot r' p{r') Aur'^dr' . 



(Al) 



The total energy of the halo within the virial radius is the 
sum, En = Wh + Th, of the potential and kinetic energies. 
The self-binding energy of the material within the virial ra- 
dius, rvir, is 



W^H(rvir 



(r') p(r') A-Kr'^ dr' 



1 / 



G 

' 2 



M{r'f , , Af^(r.i,) 
— dr H 



(A2) 



where ^(r) is the gravitational potential. Assuming hydro- 
static equilibrium with an isotropic velocity dispersion cr(r), 
the corresponding kinetic energy of material inside the virial 
radius can be expressed as 



rH(rvir) 



^cr^(r') p(r') 47rr''^ dr' . 



(A3) 



With the same assumptions, the velocity dispersion obeys 
the Jeans equation d{pa^)/dr = —pGM{r)/r^ . Provided 



2tv 



GM{r')p{r')r' dr' 



(A4) 



For our standard case of halos with the NFW density 
profile, equation (3.8), we simply integrate the Jeans equa- 
tion out to r = oo to derive u{r), assuming that the NFW 
profile an d hyd rostatic equilibrium apply at all radii (Cole 
& Lacey ( [L996| ), equation (2.14)). This is an approximation, 
since in principle we should not expect the NFW halo model 
to be valid beyond the virial radius, where material is still 
infalling. However, the velocity dispersion within the halo 
derived in this way using the NFW model has been found 
to be in good agreement with nume rical simulations (e.g. 
figures 4, 5 and 6 of Cole & Lacey 1996). Note also that 



truncating the NFW profile at the virial radius implies that 
2rH(rvir) + VKH(''vir) 7^ 0. If the integrals in equations (A2) 



and (A4) were extended to r = oo then the NFW halo model 
would exactly satisfy the virial theorem, but for the trun- 
cated model 2r(rvir)/|V[^(rvir)| is slightly greater than unity 
and varies slowly with the NFW scalelength, ankw- This 
behavi our w as also found for the N-body halos in Cole & 
Lacey (199(:) and our definitions are fully consistent with 



the way in which they defined the spin parameter Ah- 

Inserting the above definitions of JH(rvir) and EYi{r^i^) 
into equation (3.6) for Ah defines the coefficient A(aNFw) in 
the relation 



Kot = j4(aNFw)AHVH, 



(A5) 



where Kh = (GM/rvir)^''^ is the circular velocity of the halo 
at the virial radius. For the limited range 0.03 < flNFW < 0.4 
the result is well fit by yl(aNFw) ~ 4.1 -I- 1.8a^''p^. 

For the non-standard case of an isothermal density pro- 
file for the halo (with or without a core radius), we follow a 
slightly difi'erent approach. If, as above, the kinetic energy, 
Tii(r\ir), is calculated by integrating the Jeans equation to 
derive cr(r), then 2rH(rvir)/|WH(rvir)| is found to be consid- 
erably greater than unity. The discrepancy is large because 
we have extrapolated the halo density profile beyond the 
virial radius with a model whose mass does not rapidly con- 
verge. Thus, for these pro files, we prefer to define the coeffi- 
cient, A, in equation (A5) by evaluating the binding energy, 
expression (A.2), with the appropriate density profile but 
then assuming the virial theorem, 2rH(rvir)/|WH(T'vir)| = 1, 
to estimate the kinetic energy and hence the total energy 
Ey_{r Mir). This is then identical to the assumption made in 
Mo et al. (1998a) to define the energy and spin parameter. 



In the range 0.01 < a < 0.4 the resulting dependence is well 
fit by A 3.66 - 0.83 a. 



APPENDIX B: STAR FORMATION 



The set of coupled differential equations,( 4.6- 4.11), describ- 
ing an episode of star formation has the following analytic 
solutions. The mass of gas that has cooled and been accreted 
in a time t since the start of the timestep is 



AMa, 



A/cool t. 



(Bl) 
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The increase in the mass of long lived stars 
rO 1 — -R 



— McoolTcft 



1-R + P 
1~R 



exp(-t/Tcff)] 



[1 - t/Tcff - exp(-t/roff)] (B2) 



where Tcff = — -R + /3) .In terms of these quantities 

the changes in the masses of cold and hot gas are 

1- R + 13 



AMcold = AMacc 

and 

AMhot = -AMacc + 



R 



-AM* 



(B3) 



1~R 

The corresponding changes in the masses of metals are 



(B4) 



am£h = amL 



1-R 



AM* 



^"^ + ^-AMf (B5) 



R 
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R 



where 
AMfcc 



Mcool^hot t 



(B6) 
(B7) 



and 



AMf = ^ " ^ 



; (a^cZ [1 - exp(-t/rcff)] 



1 - i? + /3 ' 

-McoolTcff^hot [1 — t/TcB — CXp( — t/rcff)] 

-r=^^ [!-(!+ V^cff ) exp(-t/rcff)] 

-McooiT-cff [2-t/Tcff-(2+t/roff)exp(-t/Tcff)] 



(B8) 



For the case where there is no supply of cooling gas, 
Mcooi = 0, the above equations show that when t 2> Tcff , 
the mean metallicity of the stars that have formed is 

^* — ^cold + 



1-R + p- 

Thus our model of star formation and feedback produces an 
effective yield Pefr = (1 — e)p/(l — R + (5) which, through /3, 
and possibly also e, is a function of the potential well depth 
of the galaxy disk or bulge in which the star formation is 
occuring. 



APPENDIX C: ADIABATIC CONTRACTION 
OF HALO, DISK AND SPHEROID 

This appendix describes how we use the adiabatic contrac- 
tion model to calculate the dynamical equilibrium of the 
disk, bulge and halo. The outputs from the calculation are 
the disk and bulge radii and the halo density profile after 
deformation by the gravity of the disk and spheroid. 

CI Adiabatic Contraction of Halo 

To model the effect on the halo density profile of a galaxy 
condensing at its centre, we start by assuming that baryons 
and dark matter have the same initial density distribution, 
with total mass profile, MhoC^o) (e.g. given by the NFW pro- 
file). A fraction 1 — /h of the total mass condenses to form a 
galaxy at the centre of the halo, leaving a fraction /h of the 



mass still in the halo component. This fraction includes any 
baryons that have not yet cooled and also satellite galax- 
ies. For simplicity, any baryons left in the hot component 
are assumed to be distributed like the dark matter. We now 
assume that in response to the gravity of the disk and the 
bulge, each shell of halo matter adjusts its radius to con- 
serve its pseudo specific angular momentum rVdr), i.e. that 
rVti{r) is an adiabatic invariant. Thus, 



roKo(ro) = rVc(r) 



(CI) 



where Vdr) is the total circular velocity at radius r, ro is 
the radius of a shell before condensation of the galaxy, and 
r is the final radius of the same shell after condensation of 
the galaxy. The initial and final halo masses interior to the 
shell are related by 



Afe(r-) = /HMHo(ro), 



(C2) 



where MH(r) is the final mass halo profile. For the purpose 
of computing the circular velocity of the halo (averaged over 
spherical shells), we treat the mass distribution (including 
the disk) as spherical. This should be a better approximation 
for estimating the gravitational influence of the disk on the 
halo, than using the circular velocity due to the disk in the 
disk plane, which is somewhat larger. Thus, 



V^{r) = G [Mu{r) + Mu{r) + Msir)] /r, 



(C3) 



where Muir) and MB(r) are the disk and bulge masses in- 
terior to radius r. For consistency, the total masses should 
be related by Mho = /hA/ho + Md -I- Mb, so that t he outer 
radius of the halo is unchanged. Combining (CI), (C2) and 
(|C3|) we have. 



roMuoiro) = r [fHMuoiro) + Moir) + A/s(r)] , 



(C4) 



which relates the final radius of any halo shell to its ini- 
tial radius, once the galaxy disk and bulge profiles, Mnir) 
and MB{r), are known. The accuracy of this appro ach has 
recently been tested in Navarro & Steinmetz (2000). 



C2 Dynamical Equilibrium of Disk and Bulge 

Application of the galaxy rules described in Section give 
the total mass, Mb, and specific angular momentum, ju, 
of a galaxy disk, but, in order to use (|C4[), we require the 
complete mass profile, Alnir). To obtain this we make the 
following simplifying assumptions. First, we assume that all 
disks are well-described by an exponential surface density 
profile. 



SD(r) = exp(-r//iD), 
for which 

Mu{r) = Md [1 - (1 + r/hj,) exp{~r/hu) 



(C5) 



(C6) 



The half- mass radius, ru, is related to the scalelength, ho, 
by J^D ~ 1.68/iD. We also assume that the specific angular 
momentum of the disk is given by 



jo = kurr)VcD{rD), 



(C7) 



where Vcd(t'd) is the circular velocity in the disk plane at the 
disk half-mass radius and ko is a constant. We adopt fco = 
1.19 as is appropriate for a flat rotation curve. The value of 



© 0000 RAS, MNRAS 000, 000-000 



40 S. Cole et al. 



ko is only weakly dependent on the assumed rotation curve. 
For example, if VcD{r) is taken to be the circular velocity 
of a self-gravitating exponential disk (Binney & Tremaine 
1987| eq. 2-169), then ku = 1.09, while if V'cD(r) oc l/r^/^ 
then fcD = 1.03. 

The radius of the disk is then related to its angular 
momentum by 

jo = k^rl,Vi^j){ro) 



— kuGro 



fuMmiroo) + -fchAfo + Afe(rD) 



(C8) 



where roo is the initial radius of the shell whose final ra- 
dius is ru- The factor fch arises from the disk geometry; if 
the disk contribution to Vco (r) is computed in the spherical 
approximation, feh = 1, but here instead we calculate the 
circular velocity in the midplane of the exponential disk, us- 
ing equation (2-169) from Binney & Tremaine (1987), giving 
fch = 1.25. 

The disk half-mass radius, ru, must satisfy this equation 
and also equation (C4) evaluated at ru, i.e. 



rDoAfHo(''Do) = 



ru 



fuMaoiruo) + -Mu + AfB(rD) 



Using (OS), this can be written as 
1 



rDoMHoiroo) 



klG 



{kh - lyoMo. 



(C9) 



(CIO) 



To derive the size of the spheroidal component of a 
galaxy, we assume that the projected density profile is well 
described by t he de Vaucouleurs r^/^-law (e.g. Binney & 

, of the r^/*- 



1987, (§1-13)). The effective radius, r^ 



Tremaine 

law (the radius that contains half the mass in projection) is 
related to the half-mass radius, tb, by re = 1.35re. We de- 
fine a pseudo-specific angular momentum for the spheroid: 



jB = '■BVc(rB) 



(Cll) 



where Vc(r) is the circular velocity at r. This pseudo-specific 
angular momentum, j'b, is assumed to be conserved, except 
during galaxy mergers, when its value is determined by the 
properties of the merger remnant (see Section 4.4.2). 



This model leads to the following two equations for the 
bulge radius, in analogy to equations ( |C8[ ) and (CIO) for the 
disk radius, 



ji = rly^ir^) 



/HAfHo(rBo) + Mu{rB) + ^Mb 



and 



(C12) 



(C13) 



rBoA/Ho(»"Bo) = 77. 

Comparin g thi s last equation to the analogous equation for 
the disk, (CIO), we see that the 2nd term on the RHS has 
vanished. This results from assuming that the mean effect of 
the disk on a spherical shell of the spheroid can be estimated 
by spherically averaging the disk. 

To com pute th e disk and bulge radii we must solve equa- 
tions (ra), (|C10|), (|C1^) and (|C1^) given the disk and bulge 



masses, Mu and A/b, the initial halo profile, M^o{r\{), and 
specific angular momenta, ju and jb . The two pairs of equa- 
tions are coupled, but with care they can be solved with a 
simple iterative scheme. 



C3 Adiabatic Adjustment following Mergers or 
Disk Instabilities 

When a spheroid is formed by a galaxy merger, w e firs t 
calculate the radius of the new spheroid rnow from (4.19). 
We then compute j'b = rB.Vc{rB), with rs = rnow and 
K'(rB) = G(A^i + A^2)/2rnew. Then, using this value of 
jB , we calculate the value of tb given by the adiabatic con- 
traction model for the disk-bulge-halo equilibrium, and take 
this be the true value of tb . If a spheroid is formed via disk 
instability, we follo w the same procedure, starting from rnow 
given by equation ( 1.22 ), and assuming that the dark mat- 
ter mass involved in the initial calculation of 14^ (rs) is twice 
that within the initial half-mass radius of the galaxy. 
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